R, plyr, со сложной функцией - PullRequest
Новая
       7

R, plyr, со сложной функцией

2 голосов
/ 30 августа 2011

У меня есть следующий набор данных (CEU): 

group  x      y
1     -23     100
1     -0.90   69.62
1     -0.90   72.03
2     -23     100
2      0.69   48.01
2      0.69   45.63

Для каждого значения группы я хочу применить функции, отмеченные ниже, к каждому подмножеству значений x и y.Затем я хотел бы объединить все результаты и записать их в таблицу для экспорта.

Я не уверен, как именно применить функцию plyr для этого ... если это действительно правильный курс действий. 

x<-c(-23.0000,-0.9031,-0.9031)
y<-c(100,85.72,86.65)

par<-c(16.88,100.28,-.75,4.129)

dcrit<-function(d) { 
    sumsq<-0
    for (i in 1:length(x)){
      sumsq<-sumsq+ (y[i]-(par[1]+(par[2]-par[1])/(1+10^((x[i]-par[3])*d))))^2      
    }
    sumsq
}

S<-function(par) { 
    a<-par[1]
    b<-par[2]
    c<-par[3]
    d<-par[4]
    sumsq<-0
    for (i in 1:length(x)){
      sumsq<-sumsq+ (y[i]-(a+(b-a)/(1+10^((x[i]-c)*d))))^2      
    }
    sumsq
}
optim(par,S)

CEU <- read.csv(file="C:/files/CEU.csv",head=TRUE,sep=",")
CEU

data <- ddply(CEU,.(group),function(xy) 
{
par[1]<-min(y)
par[2]<-100
par[3]<-x[[which.min(abs(y-50))]]
par[4]<-optimize(dcrit,interval=c(-100,100))$minimum

o<-optim(par,S)
par<-o$par

a<-par[1];
b<-par[2];
c<-par[3];
d<-par[4];

k<-(b-a)/(20-a)-1
if (k>0) ec20<-c+1/d*log10(k) else ec20<-NA
ec20

z<-(b-a)/(50-a)-1
 if (z>0) ec50<-c+1/d*log10(z) else ec50<-NA
ec50

j<-(b-a)/(80-a)-1
if (j>0) ec80<-c+1/d*log10(j) else ec80<-NA
ec80

data.frame(ec20, ec50, ec80)

})

data

Код выполняется без ошибок, но только наисходные значения x и y, установленные с помощью: 

 x<-c(-23.0000,-0.9031,-0.9031)
 y<-c(100,85.72,86.65)

Значения x и y в наборе данных CEU не используются ddply.Они не заменяют исходные x и y итеративным способом, как это делают значения групп.данные имеют соответствующее количество групп, и значения ec20 / ec50 / ec80 являются точными, но только для исходных x и y.

> data
   group       ec20       ec50       ec80
1      1 -0.3652977 -0.6843279 -0.8530892
2      2 -0.3652977 -0.6843279 -0.8530892
3      3 -0.3652977 -0.6843279 -0.8530892
4      4 -0.3652977 -0.6843279 -0.8530892
5      5 -0.3652977 -0.6843279 -0.8530892

1 Ответ

3 голосов
/ 30 августа 2011

Похоже, у вас все правильно, вам просто нужно произвести вывод.

Я предполагаю, что это где ваши выходы? 

k<-(b-a)/(20-a)-1
if (k>0) ec20<-c+1/d*log10(k) else ec20<-NA
ec20

z<-(b-a)/(50-a)-1
 if (z>0) ec50<-c+1/d*log10(z) else ec50<-NA
ec50

j<-(b-a)/(80-a)-1
if (j>0) ec80<-c+1/d*log10(j) else ec80<-NA
ec80

Положите их в data.frame в конце функции: 

    ...
    data.frame(ec20, ec50, ec80)
}

Теперь вы получите data.frame со всеми из них, с тремя столбцами для ec20, ec50 и ec80

К вашей проблеме с optim: я думаю, что проблема заключается в 

par[3]<-x[which.min(abs(y-50))]

Single [ в R не является обычным индексом - он получает фрагмент - в этом случаеdata.frame столбцы.Эта линия превращает par из числового вектора в list.Добавьте больше скобок: 

par[3]<-x[[which.min(abs(y-50))]]
...