Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
################クレームモデル###########################
##各分布の密度関数
layout(matrix(1:6,nrow=2,byrow=FALSE))
plot(dpois(0:20,5))
plot(dnbinom(0:20,5,1/2))
curve(dgamma(x,2,0.01),0,500)
curve(dnorm(x,200,50),0,500)
curve(dlnorm(x,4.7,1.094),0,500)
##クレーム総額を作る関数
SCC=function(n,type1,para11,para12,type2,para21,para22){
#頻度はポアソン(1)または負の二項分布(幾何分布含む)(2)に従う。
#クレーム額はガンマ分布(指数分布含む)(1)、正規分布(2)、対数正規分布(3)に従う
result=rep(0,n) #長さnのベクトルに各試行のクレーム総額を入れていく
for(i in 1:n){
#クレーム件数を作る
if(type1==1){
N=rpois(1,para11) #ポアソン分布の乱数
}else{
N=rnbinom(1,para11,para12) #負の二項分布の乱数
}
#クレーム額を作る
if(type2==1){
cc=rgamma(N,para21,para22)
}else if(type2==2){
cc=rnorm(N,para21,para22)
cc[cc<0]=0
}else{
cc=rlnorm(N,para21,para22)
}
result[i]=sum(cc)
}
##経験分布をヒストグラムで書かせる
hist(result,freq = FALSE, breaks = 20,right=FALSE)
##標本平均を赤線で引く
abline(v=mean(result),col=2)
##1σ分ずらして緑線を引く
abline(v=c(mean(result)+(var(result))^(1/2),mean(result)-(var(result))^(1/2)),col=3)
}
##クレーム総額サンプル
windows()
layout(matrix(1:6,nrow=2,byrow=FALSE))
SCC(1000,1,5,0,1,2,0.01)         #ポアソン分布とガンマ分布
SCC(1000,2,5,1/2,1,2,0.01) #負の二項分布とガンマ分布
SCC(1000,1,5,0,2,200,50) #ポアソン分布と正規分布
SCC(1000,2,5,1/2,2,200,50) #負の二項分布と正規分布
SCC(1000,1,5,0,3,4.7,1.094) #ポアソン分布と対数正規分布
SCC(1000,2,5,1/2,3,4.7,1.094) #負の二項分布と対数正規分布
##NB(1,0.1)とGamma(1,0.1)の複合分布を考える
windows()
SCC(1000,2,1,0.1,1,1,0.1)     #経験分布を書く
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment