損保数理 日本アクチュアリー会

#複合ポアソン分布　2-15

library(dplyr)

n=10

x=sample(1:100,n,replace=T)

p=x/sum(x);lam=sum(x)

t=sample(1:100,n,replace=T)

(lam*p*(exp(t)-1))

lam_i=t*lam*p

V_S=lam*sum(p*t^2);E_S=sum(t*p);

under=E_S+qnorm(0.05/2)*sqrt(V_S)

upper=E_S+qnorm(1-0.05/2)*sqrt(V_S)

#Jung method 4-13

I=5;J=3

nij=array(sample(1:10,I*J,replace=T),dim=c(I,J))

rij=array(sample(1:10,I*J,replace=T),dim=c(I,J))

#Broyden Quasi-Newton Algorithm p.213

f=function(s){

x=s[1:I];y=s[-c(1:I)]

f1=apply(-nij*c(y)+t(t(nij*rij)/c(x)),1,sum)

f2=apply(-t(t(nij)*c(x))+((nij*rij)/c(y)),2,sum)

return(c(f1,f2))

}

X=rep(1,I+J)

ite=1000

H=diag(f(X))

eta=0.001

for(l in 1:ite){

if(sum(abs(f(X)))>10^(-9)){  

X_pre=X  

X=X-eta*H%*%f(X)  

s=X-X_pre

y=f(X)-f(X_pre)

H=H+((s-H%*%y)/as.numeric(t(s)%*%H%*%y))%*%t(s)%*%H

print(sum(abs(f(X))))

}

}