1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
|
interpweights <- function(w, v1, v2){
#Given L=(w,v1), compute neww such that newL=(new,v2)=L in distribution
cumw <- cumsum(w)
neww <- splinefun(v1, cumw, method= "monoH.FC")(v2, deriv=1)
#neww <- diff(newcumw)
interpweights <- neww/sum(neww)
return(interpweights)
}
interpvalues <- function(w, v, neww){
## Given a distribution D=(w,v), compute new values
## such that Dnew=(neww, newv) equals D in distribution
cumw <- cumsum(w)
cdf <- splinefun(v, cumw, method="hyman")
eps <- 1e-3
newv <- rep(0, length(neww))
cumneww <- cumsum(neww)
mid <- 0
for(i in 1:length(neww)){
iter <- 0
## do binary search
hi <- cdf(1)
lo <- mid
if(hi < cumneww[i]){
newv[i] <- hi
next
}
if(cdf(lo) > cumneww[i]){
newv[i] <- lo
next
}
mid <- (lo+hi)/2
iter <- 0
while(abs(cdf(mid) - cumneww[i])>eps){
if(cdf(mid) > cumneww[i]){
hi <- mid
}else{
lo <- mid
}
mid <- (lo+hi)/2
}
newv[i] <- mid
}
return(newv)
}
interpvalues.distr <- function(w, v, neww){
## same as interpvalues, but using the distr
## package. need to check how good it is
require(distr)
D <- DiscreteDistribution(v, w)
return( q(D)(cumsum(neww)) )
}
adjust_scenario <- function(scenario, epsilon){
1-(1-scenario)^(1/(1+epsilon))
}
adjust_weights <- function(weights, scenario, epsilon){
interpweights(weights,scenario,adjust_scenario(scenario,epsilon))
}
obj <- function(epsilon, vecpv, prob, support, cte){
newprob <- adjust_weights(prob, support, epsilon)
return( 1 - crossprod(newprob, vecpv) - cte)
}
tweak <- function(min, max, vecpv, prob, support, cte){
mid <- (min + max)/2
objective <- obj(mid, vecpv, prob, support, cte)
while( abs(objective) > 1e-6){
if(objective > 0){
min <- mid
}else{
max <- mid
}
mid <- (min+max)/2
objective <- obj(mid, vecpv, prob, support, cte)
}
return( mid )
}
interpweightsadjust <- function(w, v1, v2, vecpv){
interpweightsadjust <- interpweights(w, v1, v2)
epsilon <- tweak(-0.5, 0.5, vecpv, interpweightsadjust, v2, 1)
return( adjust_weights(interpweightsadjust, v2, epsilon) )
}
transformweightslike <- function(p1, v1, p2, v2, p, v){
cump2 <- cumsum(p2)
cump1 <- cumsum(p1)
P1 <- splinefun(v1,cump1,method= "monoH.FC")
dP1 <- function(x){P1(x,deriv=1)}
pomme <- interpweights(p2,v2,v)
pomme <- cumsum(pomme)
r <- rep(0,length(pomme))
for(i in 1:length(pomme)){
r[i] <- inverse(P1,dP1,pomme[i])
}
return(r)
}
clipw <- function(x){
write(x,file="clipboard",sep="\n")
}
clipr <- function(){
scan(file="clipboard")
}
sclipr <- function(){
scan(file="clipboard",what="character")
}
inverse <- function(f,Df,x, x0=x){
#inverse a function by the newton's method.
x1 <- x0-f(x0)/(Df(x0)-1)
counter <- 0
while(abs(x1-x0)>1e-6&&counter<500){
x0 <- x1
x1 <- x0-(f(x0)-x)/Df(x0)
counter <- counter+1
}
return(x1)
}
|
