R for LC+cohort

#导入美国数据
USA.mortality<-read.csv("/Users/liuqing/Downloads/死亡率数据库/美国/美国死亡率1950-2017.csv",header=TRUE)
USA.population<-read.csv("/Users/liuqing/Downloads/死亡率数据库/美国/美国人数1950-2017.csv",header=TRUE)
USA.deathpopulation<-read.csv("/Users/liuqing/Downloads/死亡率数据库/美国/美国死亡人数1950-2017.csv",header=TRUE)
USA.gdp<-read.csv("/Users/liuqing/Downloads/死亡率数据库/美国/女性/美国人均gdp.csv",head=TRUE)#1960-2018
m_all<-matrix(rep(0,100*68),100,68)   #死亡率
lnm_all<-matrix(rep(0,100*68),100,68) #ln死亡率
D_all<-matrix(rep(0,100*68),100,68)   #计算的死亡人数
E_all<-matrix(rep(0,100*68),100,68)   #总人数
d_all<-matrix(rep(0,100*68),100,68)   #实际观测的死亡人数
#女性1950-2017年的所有数据
for (t in 1:68){for(x in 1:100){E_all[x,t]<-USA.population[x+100*(t-1),3];   #总人数，100行，68列m_all[x,t]<-USA.mortality[x+100*(t-1),3];    #女性死亡率D_all[x,t]<-USA.deathpopulation[x+100*(t-1),3] #实际观测的死亡人数lnm_all[x,t]<-log(m_all[x,t]);}
}
#用于实验的数据
E<-E_all[,1:63] #1950-2012
m<-m_all[,1:63] #1950-2012
lnm<-lnm_all[,1:63]
m_original<-lnm_all[,64:68]#用于预测的死亡率
y<-D_all[,1:63]#实际观测到的死亡人数
#添加age-period-cohort
a<-c(rep(0,100))
b0<-c(rep(1,100))
b1<-c(rep(1,100))
k<-c(rep(0,63))
l<-c(rep(0,162))
#1.1 估计的死亡率
mm<-function(a,b0,b1,k,l){m_hat<-matrix(rep(0,100*63),100,63)for(t in 1950:2012){for(x in 0:99){for(z in 1851:2012){if(z==t-x){m_hat[x+1,t-1949]<-exp(a[x+1]+b0[x+1]*k[t-1949]+b1[x+1]*l[t-x-1850])}}}}return(m_hat)
}
#1.2 估计的死亡人数
yy<-function(a,b0,b1,k,l){y_hat<-matrix(rep(0,100*63),100,63)E<-E_all[,1:63]for(x in 1:100){for(t in 1:63){y_hat[x,t]<-E[x,t]*m_hat[x,t]      #估计的死亡人数}}return(y_hat)
}#1.3 似然函数的值
RSS<-function(a,b0,b1,k,l){L<-matrix(rep(NA,100*63),100,63)y_hat<-yy(a,b0,b1,k,l)m_hat<-mm(a,b0,b1,k,l)for(x in 1:100){for(t in 1:63){L[x,t]<-2*(y[x,t]*log(y[x,t]/y_hat[x,t])-(y[x,t]-y_hat[x,t]))}}suml<-sum(L)return(suml)
}#2.1 alpha
aa<-function(a,b0,b1,k,l){y_hat<-yy(a,b0,b1,k,l)for(x in 1:100){a[x]<-a[x]+sum(y[x,]-y_hat[x,])/(sum(y_hat[x,]))}return(a)
}#2.2 kt
kk<-function(a,b0,b1,k,l){y_hat<-yy(a,b0,b1,k,l)for(t in 1:63){k[t]<-k[t]+sum((y[,t]-y_hat[,t])*b1)/(sum(y_hat[,t]*(b1^2)))}return(k)
}#2.3 beta(1)
bb1<-function(a,b0,b1,k,l){y_hat<-yy(a,b0,b1,k,l)for(x in 1:100){b1[x]<-b1[x]+sum((y[x,]-y_hat[x,])*k)/(sum(y_hat[x,]*(k^2)))}return(b1)
}#2.4 cohort
ll<-function(a,b0,b1,k,l){y_hat<-yy(a,b0,b1,k,l) for(x in 0:99){           #最后对x循环for(z in 1851:2012){    #再对z循环for(t in 1950:2012){  #先对t循环if(t-x==z){l[z-1850]<-l[z-1850]+sum((y[x+1,t-1949]-y_hat[x+1,t-1949])*b0[x+1])/(sum(y_hat[x+1,t-1949]*(b0[x+1]^2)))}}}}return(l)
}#2.5
bb0<-function(a,b0,b1,k,l){y_hat<-yy(a,b0,b1,k,l)for( x in 0:99){for(t in 1950:2012){for(z in 1851:2012){if(z==t-x){b0[x+1]<-b0[x+1]+sum((y[x+1,t-1949]-y_hat[x+1,t-1949])*l[t-x-1850])/sum(y_hat[x+1,t-1949]*(l[t-x-1850])^2)}}return(b0)}}
}
#3
j=1
diff_L<-RSS(a,b0,b1,k,l)
while(abs(diff_L)>10^(-10))
{logL0<-RSS(a,b0,b1,k,l)if(j%%5==1){a<-aa(a,b0,b1,k,l)}if(j%%5==2){k<-kk(a,b0,b1,k,l)}if(j%%5==3){b0<-bb0(a,b0,b1,k,l)}if(j%%5==4){b1<-bb1(a,b0,b1,k,l)}if(j%%5==0){l<-ll(a,b0,b1,k,l)}logL1<-RSS(a,b0,b1,k,l)diff_L<-logL1-logL0j=j+1
}

R for LC+cohort相关推荐

hdu-3071 Gcd Lcm game---质因数分解+状态压缩+线段树
题目链接: http://acm.hdu.edu.cn/showproblem.php?pid=3071 题目大意: 给定一个长度为n的序列m次操作,操作的种类一共有三种查询 L :查询一个区间的所 ...
BZOJ1251: 序列终结者
[传送门:BZOJ1251] 简要题意: 给出一个长度为n的序列,有m个操作,3种操作: 1 l r k将l到r的数增加k 2 l r将l到r的数列翻转 3 l r求出l到r的最大值题解: 裸SPL ...
【GDKOI2016Day1T1-魔卡少女】【拆位】线段树维护区间内所有连续子区间的异或和...
题意:给出N个数,M个操作.操作有修改和询问两种,每次修改将一个数改成另一个数,每次询问一个区间的所有连续子区间的异或和.n,m<=100000,ai<=1000 题解: 当年(其实也就是 ...
【BZOJ】3339: Rmq Problem 3585: mex（线段树+特殊的技巧）
http://www.lydsy.com/JudgeOnline/problem.php?id=3585 好神的题. 但是!!!!!!!!!!!!!!我线段树现在要开8倍空间才能过!!!!!!!!!! ...
bzoj2243 [SDOI2011]染色
树链剖分,线段树维护三个值:区间颜色数.区间左端颜色.区间右端颜色,当左区间右端颜色=右区间左端颜色时候,总颜色-1,树上的维护类似. 代码 1 #include<cstdio> 2 co ...
BZOJ 1014 [JSOI2008]火星人prefix
splay维护hash值: 看到大佬们都在打数据结构,我说好的不打数据结构又自己打脸了. 为了写这个昨天还特意去打了Splay的普通平衡树,自从我学会Treap以来第一次用splayA掉普通平衡树QA ...
BZOJ4373: 算术天才⑨与等差数列
Description 算术天才⑨非常喜欢和等差数列玩耍. 有一天,他给了你一个长度为n的序列,其中第i个数为a[i]. 他想考考你,每次他会给出询问l,r,k,问区间[l,r]内的数从小到大排序后能 ...
BZOJ 2333 【SCOI2011】棘手的操作
题目链接:棘手的操作网上的题解大部分都是在线用可并堆艹--但是树高严格\(\log\)的可并堆我不会啊--还是离线大法好-- 我们可以先把所有的合并操作用并查集给处理好,把得到的森林记录下来.然后, ...
Codeforces Round #395 (Div. 2)（未完）
2.2.2017 9:35~11:35 A - Taymyr is calling you 直接模拟 #include <iostream> #include <cstdio> ...

R for LC+cohort

R for LC+cohort

R for LC+cohort相关推荐

最新文章

热门文章