Версия для печати темы
Форум врачей-аспирантов _ Медицинская статистика _ Подскажите со стат.методом
Автор: Де бин Анатолий 12.10.2017 - 13:33
Приветствую. , пожалуйста, мне тут метод выбрать в такой задаче.
У меня есть данные по 80 приборам, по каждому циклу, для каждого прибор, разные датчики дают указание. Как я могу предложить на основании этих данных, что прибор неисправен? Какой статистический метод я должен использовать?
Прибор считается вышедшим из строя, если больше измерений по нему нет. Например, есть измерения по циклам 1-151. Значит, что на 152 цикле прибор вышел из строя.
Может ли тут помочь анализ выживаемости?
Автор: leo_biostat 14.10.2017 - 19:24
Анатолий, !
У Вас отличный массив данных! 25 признаков и 16139 наблюдений. Задача непростая, и поэтому не стоит ориентироваться на элементарные, примитивные методы. Лет 40 назад у меня была одна подобная задача, когда я работал в одном техническом НИИ.
Нужно было сравнить полтора десятка приборов, и установить идентичность и различие этих приборов. Число наблюдений было гораздо больше. И тоже возникла проблема выбора оптимальной коллекции методов анализа. Исходя из своего опыта решения этой и иных подобных задач, рекомендую использовать следующие методы анализа:
1. Кластерный анализ.
2. Анализ таблиц сопряжённости, с анализом структуры связи.
3. Дискриминантный анализ.
4. Анализ канонической корреляции, с графикой распределения по новым осям.
5. Анализ функций распределения признаков по имеющимся группировкам, и установленным кластерам.
6. Сравнение параметров распределения уже имеющихся группировок, и установленных кластеров.
По результатам этих анализов нужно будет установить те приборы, которые статистически значимо отличаются по своим показаниям от идентичных группировок, а также установить наличие аномальных значений измеряемых величин у конкретных приборов в конкретных группах. Это и будет аргументацией для принятия подобных приборов аномальными.
Желаю успеха!
Цитата(Де бин Анатолий @ 12.10.2017 - 13:33)
Приветствую.
, пожалуйста, мне тут метод выбрать в такой задаче.
У меня есть данные по 80 приборам, по каждому циклу, для каждого прибор, разные датчики дают указание. Как я могу предложить на основании этих данных, что прибор неисправен? Какой статистический метод я должен использовать?
Прибор считается вышедшим из строя, если больше измерений по нему нет. Например, есть измерения по циклам 1-151. Значит, что на 152 цикле прибор вышел из строя.
Может ли тут помочь анализ выживаемости?
Автор: p2004r 15.10.2017 - 20:33
Цитата(Де бин Анатолий @ 12.10.2017 - 13:33)
Приветствую.
, пожалуйста, мне тут метод выбрать в такой задаче.
У меня есть данные по 80 приборам, по каждому циклу, для каждого прибор, разные датчики дают указание. Как я могу предложить на основании этих данных, что прибор неисправен? Какой статистический метод я должен использовать?
Прибор считается вышедшим из строя, если больше измерений по нему нет. Например, есть измерения по циклам 1-151. Значит, что на 152 цикле прибор вышел из строя.
Может ли тут помочь анализ выживаемости?
Можно использовать анализ выживаемости, но тогда придется "ручками" сводить многомерную последовательность индивидуальных показаний прибора к строке фиксированной длинны.
Это задача в целом стандартная и описывается паттерном "many to one". Решение тоже стандартное -- грузиться library(keras) и описывается сетка с LSTM (хот для такого размера выборки GRU) слоем который и сворачивает историю прибора в вектор стандартной длины. На выходе предиктор "авария есть -- аварии нет". На входе список с историями приборов + всякая его нормализация.
Если "прогноз момента отказа" надо сделать, то тогда придется делать архитектуру seq2seq, и на выход модели пристраивать еще один слой GRU который будет выдавать цепочку непрерывную "все ок"кончающуюся "всё плохо" для сломавшихся приборов.
В общем тут надо конкретно разбираться что предсказывать надо.
PS
Но конечно можно и марковские цепочки попробовать пообучать, хоть это и весьма геморройное занятие.
Автор: Де бин Анатолий 16.10.2017 - 12:41
Код
Это задача в целом стандартная и описывается паттерном "many to one". Решение тоже стандартное -- грузиться library(keras) и описывается сетка с LSTM (хот для такого размера выборки GRU) слоем который и сворачивает историю прибора в вектор стандартной длины. На выходе предиктор "авария есть -- аварии нет". На входе список с историями приборов + всякая его нормализация.
установил библиотечку keras
нашел в хелпе описание
Код
layer_conv_lstm_2d {keras} R Documentation
Convolutional LSTM.
Description
It is similar to an LSTM layer, but the input transformations and recurrent transformations are both convolutional.
Usage
layer_conv_lstm_2d(object, filters, kernel_size, strides = c(1L, 1L),
padding = "valid", data_format = NULL, dilation_rate = c(1L, 1L),
activation = "tanh", recurrent_activation = "hard_sigmoid",
use_bias = TRUE, kernel_initializer = "glorot_uniform",
recurrent_initializer = "orthogonal", bias_initializer = "zeros",
unit_forget_bias = TRUE, kernel_regularizer = NULL,
recurrent_regularizer = NULL, bias_regularizer = NULL,
activity_regularizer = NULL, kernel_constraint = NULL,
recurrent_constraint = NULL, bias_constraint = NULL,
return_sequences = FALSE, go_backwards = FALSE, stateful = FALSE,
dropout = 0, recurrent_dropout = 0, batch_size = NULL, name = NULL,
trainable = NULL, weights = NULL, input_shape = NULL)
А прикладного примера в хелпе нет, обычно он есть. Мне без него сложно перенести идею на мои данные.
в роле object
мой датасет должен быть?
setwd(mydir)
mydat=read.csv(mydat)
Автор: p2004r 16.10.2017 - 19:54
Цитата(Де бин Анатолий @ 16.10.2017 - 12:41)
А прикладного примера в хелпе нет, обычно он есть. Мне без него сложно перенести идею на мои данные.
У самого кераса точно также с документацией, все лежит в руководствах.
На странице пакета в cran много руководств лежит (и в этих руководствах много ссылок на дополнительную информацию с сайта Революшен)
https://cran.r-project.org/web/packages/keras/index.html
Начинаете с Getting Started with Keras и потом Guide to the Functional API.
Автор: p2004r 18.10.2017 - 16:03
Цитата(Де бин Анатолий @ 16.10.2017 - 12:41)
Код
Это задача в целом стандартная и описывается паттерном "many to one". Решение тоже стандартное -- грузиться library(keras) и описывается сетка с LSTM (хот для такого размера выборки GRU) слоем который и сворачивает историю прибора в вектор стандартной длины. На выходе предиктор "авария есть -- аварии нет". На входе список с историями приборов + всякая его нормализация.
установил библиотечку keras
нашел в хелпе описание
Код
layer_conv_lstm_2d {keras} R Documentation
Convolutional LSTM.
Description
It is similar to an LSTM layer, but the input transformations and recurrent transformations are both convolutional.
Usage
layer_conv_lstm_2d(object, filters, kernel_size, strides = c(1L, 1L),
padding = "valid", data_format = NULL, dilation_rate = c(1L, 1L),
activation = "tanh", recurrent_activation = "hard_sigmoid",
use_bias = TRUE, kernel_initializer = "glorot_uniform",
recurrent_initializer = "orthogonal", bias_initializer = "zeros",
unit_forget_bias = TRUE, kernel_regularizer = NULL,
recurrent_regularizer = NULL, bias_regularizer = NULL,
activity_regularizer = NULL, kernel_constraint = NULL,
recurrent_constraint = NULL, bias_constraint = NULL,
return_sequences = FALSE, go_backwards = FALSE, stateful = FALSE,
dropout = 0, recurrent_dropout = 0, batch_size = NULL, name = NULL,
trainable = NULL, weights = NULL, input_shape = NULL)
А прикладного примера в хелпе нет, обычно он есть. Мне без него сложно перенести идею на мои данные.
в роле object
мой датасет должен быть?
setwd(mydir)
mydat=read.csv(mydat)
Массив входных данных вот как то так готовиться
Код
df.data <- read.csv2("mydat.csv")
library(dplyr)
library(keras)
data.list <- lapply(split(df.data, df.data$id), as.list)
pribor_order <- names(data.list)
names(data.list) <- NULL
data.list <- lapply(data.list, function(l) {names(l) <- NULL; l[5:25]})
data.list <- lapply(data.list, function(l) lapply(l, as.numeric))
data.list <- lapply(data.list, function(l) as.matrix(data.frame(l)))
В результате мы можем подготовить входные истории для GRU
Код
> str(pad_sequences(data.list))
int [1:80, 1:362, 1:21] 0 0 0 0 0 0 0 0 0 0 ...
максимальная история 362 события, остальные добиты "0" до той же длинны
Автор: p2004r 18.10.2017 - 18:22
Вот собственно сетка целиком.
Код
df.data <- read.csv2("mydat.csv")
library(dplyr)
library(keras)
data.list <- lapply(split(df.data, df.data$id), as.list)
pribor_order <- names(data.list)
names(data.list) <- NULL
data.list <- lapply(data.list, function(l) {names(l) <- NULL; l[5:25]})
data.list <- lapply(data.list, function(l) lapply(l, as.numeric))
data.list <- lapply(data.list, function(l) as.matrix(data.frame(l)))
pribor.out <- sapply(pribor_order, function(n) sum(df.data$id == n))
seqLength <- 50
in_seq_event <- pad_sequences(data.list, truncating = 'post', maxlen = seqLength)
## model
drop <- 0.4
K <- backend()
K$clear_session()
main_input <- layer_input(shape = c(seqLength, 21), dtype = 'float32', name = 'main_input')
data_in_norm <- main_input %>%
layer_batch_normalization()
## bidirectional()
lstm_out.forw <- data_in_norm %>%
layer_lstm(units = 8,
dropout = 0.2,
recurrent_dropout = 0.2) %>%
layer_batch_normalization()
lstm_out.back <- data_in_norm %>%
layer_lstm(units = 8,
dropout = 0.2,
recurrent_dropout = 0.2,
go_backwards=T) %>%
layer_batch_normalization()
main_output <- layer_concatenate(c(lstm_out.forw,
lstm_out.back)) %>%
layer_dense(units = 8, activation = 'relu') %>%
layer_batch_normalization() %>%
layer_dropout(drop) %>%
layer_dense(units = 8, activation = 'relu') %>%
layer_batch_normalization() %>%
layer_dropout(drop) %>%
layer_dense(units = 1, activation = 'linear', name = 'main_output')
model <- keras_model(
inputs = c(main_input),
outputs = c(main_output)
)
model %>%
compile(
optimizer = 'adam',
loss = 'mean_squared_error'
)
summary(model)
history <- model %>%
fit(x = list(in_seq_event),
y = list(pribor.out),
validation_split = 0.1,
epochs = 340,
batch_size = 8)
plot(history)
pribor_predict <- model %>%
predict(x = list(in_seq_event),
batch_size = 8,
verbose = 1)
plot(pribor.out, pribor_predict)
По первым 50 записям логов приборов
Автор: p2004r 19.10.2017 - 10:22
Нормализацию забыл, вот так намного устойчивее и точнее...
В принципе приборов похоже две группы.
Код
df.data <- read.csv2("mydat.csv")
str(df.data)
df.data.norm <- (data.frame(lapply(df.data[,!sapply(df.data, is.factor)], scale))+2)[, -c(1,2,3,4,10,12,13)]
library(dplyr)
library(keras)
data.list <- lapply(split(df.data.norm, df.data$id), as.list)
pribor_order <- names(data.list)
names(data.list) <- NULL
data.list <- lapply(data.list, function(l) {names(l) <- NULL; l})
data.list <- lapply(data.list, function(l) lapply(l, as.numeric))
data.list <- lapply(data.list, function(l) as.matrix(data.frame(l)))
## data.list <- lapply(data.list, normalize)
pribor.out <- sapply(pribor_order, function(n) sum(df.data$id == n))
seqLength <- 50
in_seq_event <- pad_sequences(data.list, truncating = 'post', maxlen = seqLength)
## model
drop <- 0.4
K <- backend()
K$clear_session()
main_input <- layer_input(shape = c(seqLength, 18), dtype = 'float32', name = 'main_input')
data_in_norm <- main_input %>%
layer_batch_normalization()
## bidirectional()
lstm_out.forw <- data_in_norm %>%
layer_gru(units = 8,
dropout = 0.2,
recurrent_dropout = 0.2,
return_sequences = T) %>%
layer_batch_normalization() %>%
layer_gru(units = 8,
dropout = 0.2,
recurrent_dropout = 0.2) %>%
layer_batch_normalization()
lstm_out.back <- data_in_norm %>%
layer_gru(units = 8,
dropout = 0.2,
recurrent_dropout = 0.2,
return_sequences = T,
go_backwards=T) %>%
layer_batch_normalization() %>%
layer_gru(units = 8,
dropout = 0.2,
recurrent_dropout = 0.2,
go_backwards=T) %>%
layer_batch_normalization()
main_output <- layer_concatenate(c(lstm_out.forw,
lstm_out.back)) %>%
layer_dense(units = 8, activation = 'relu') %>%
layer_batch_normalization() %>%
layer_dropout(drop) %>%
# layer_dense(units = 8, activation = 'relu') %>%
# layer_batch_normalization() %>%
# layer_dropout(drop) %>%
layer_dense(units = 8, activation = 'relu') %>%
layer_batch_normalization() %>%
layer_dropout(drop) %>%
layer_dense(units = 1, activation = 'linear', name = 'main_output')
model <- keras_model(
inputs = c(main_input),
outputs = c(main_output)
)
model %>%
compile(
## optimizer = 'rmsprop',
## optimizer = 'nadam',
optimizer = 'adam',
loss = 'mean_squared_error'
## loss = 'binary_crossentropy'
## loss_weights = c(1.0, 0.2)
)
summary(model)
history <- model %>%
fit(x = list(in_seq_event),
y = list(pribor.out),
## y = list(labels, labels),
validation_split = 0.1,
#class_weight = list("0"=0.1392319, "1"=1),
epochs = 177*2+77,
#callbacks = list(callback_reduce_lr_on_plateau(monitor = "val_loss", factor = 0.1),
# callback_tensorboard("logs/run_a", write_images = T, embeddings_freq = 1)),
batch_size = 8)
plot(history)
pribor_predict <- model %>%
predict(x = list(in_seq_event),
batch_size = 8,
verbose = 1)
plot(pribor.out, pribor_predict)
Автор: p2004r 21.10.2017 - 16:03
Ну и собственно кластеризация по "импульсным переходным функция" восстановленным рекуррентными слоями модели
Код
lstm_out.forw_model <- keras_model(inputs = model$input,
outputs = lstm_out.forw)
lstm_out.forw_output <- predict(lstm_out.forw_model, x = list(in_seq_event), batch_size = 8, verbose = 1)
lstm_out.back_model <- keras_model(inputs = model$input,
outputs = lstm_out.back)
lstm_out.back_output <- predict(lstm_out.back_model, x = list(in_seq_event), batch_size = 8, verbose = 1)
str(lstm_out.forw_output)
str(lstm_out.back_output)
pairs(cbind(prcomp(cbind(lstm_out.forw_output,
lstm_out.back_output),
center=T, scale.=T)$x[,1:3],
pribor.out))
Как видно из графика действительно две группы приборов + некая функция "наработки на отказ".
Автор: p2004r 21.10.2017 - 16:08
Ну и по именам приборы в разных группах
Код
> pribor_order[which(prcomp(cbind(lstm_out.forw_output,
+ lstm_out.back_output),
+ center=T, scale.=T)$x[,1] < 0)]
[1] "Engine_1" "Engine_10" "Engine_11" "Engine_15" "Engine_17" "Engine_18"
[7] "Engine_2" "Engine_28" "Engine_3" "Engine_30" "Engine_31" "Engine_33"
[13] "Engine_41" "Engine_43" "Engine_44" "Engine_45" "Engine_46" "Engine_47"
[19] "Engine_48" "Engine_5" "Engine_51" "Engine_54" "Engine_58" "Engine_59"
[25] "Engine_61" "Engine_64" "Engine_65" "Engine_67" "Engine_69" "Engine_7"
[31] "Engine_72" "Engine_74" "Engine_78" "Engine_79" "Engine_9"
>
> pribor_order[which(prcomp(cbind(lstm_out.forw_output,
+ lstm_out.back_output),
+ center=T, scale.=T)$x[,1] >= 0)]
[1] "Engine_12" "Engine_13" "Engine_14" "Engine_16" "Engine_19" "Engine_20"
[7] "Engine_21" "Engine_22" "Engine_23" "Engine_24" "Engine_25" "Engine_26"
[13] "Engine_27" "Engine_29" "Engine_32" "Engine_34" "Engine_35" "Engine_36"
[19] "Engine_37" "Engine_38" "Engine_39" "Engine_4" "Engine_40" "Engine_42"
[25] "Engine_49" "Engine_50" "Engine_52" "Engine_53" "Engine_55" "Engine_56"
[31] "Engine_57" "Engine_6" "Engine_60" "Engine_62" "Engine_63" "Engine_66"
[37] "Engine_68" "Engine_70" "Engine_71" "Engine_73" "Engine_75" "Engine_76"
[43] "Engine_77" "Engine_8" "Engine_80"
Форум Invision Power Board (http://www.invisionboard.com)
© Invision Power Services (http://www.invisionpower.com)