Diff of /src/modules/glm/samplers/ConjugateFMethod.cc [01a998] .. [081fbe]  Maximize  Restore

Switch to unified view

a/src/modules/glm/samplers/ConjugateFMethod.cc b/src/modules/glm/samplers/ConjugateFMethod.cc
...
...
70
        coef[i] += getScale(schildren[i], chain);
70
        coef[i] += getScale(schildren[i], chain);
71
    }
71
    }
72
    gv->setValue(&xold, 1, chain);
72
    gv->setValue(&xold, 1, chain);
73
}
73
}
74
74
75
namespace glm {
75
76
76
ConjugateFMethod::ConjugateFMethod(GraphView *gv1, GraphView *gv2, 
77
    ConjugateFMethod::ConjugateFMethod(GraphView *gv1, GraphView *gv2, 
77
                   unsigned int chain)
78
                       unsigned int chain)
78
    : _gv1(gv1), _gv2(gv2), _chain(chain),
79
  : _gv1(gv1), _gv2(gv2), _chain(chain),
79
      _scale(1), _tau(gv1->nodes()[0]->value(chain)[0]),
80
    _scale(1), _tau(gv1->nodes()[0]->value(chain)[0]),
80
      _coef(0)
81
    _coef(0)
81
{
82
    if(!_gv1->deterministicChildren().empty() && checkScale(_gv1, true)) 
83
    {
82
    {
83
  if(!_gv1->deterministicChildren().empty() && checkScale(_gv1, true)) 
84
  {
84
    //One-off calculation of fixed scale transformation
85
        //One-off calculation of fixed scale transformation
85
    _coef = new double[_gv1->stochasticChildren().size()];
86
        _coef = new double[_gv1->stochasticChildren().size()];
86
    calCoef(_coef, _gv1, chain);
87
        calCoef(_coef, _gv1, chain);
88
  }
89
    }
87
    }
90
    
88
}
89
90
ConjugateFMethod::~ConjugateFMethod()
91
    ConjugateFMethod::~ConjugateFMethod()
91
{
92
    {
92
    delete [] _coef;
93
  delete [] _coef;
93
}
94
    }
94
95
    
95
96
    
96
//#include <iostream> //debuggin
97
    //#include <iostream> //debuggin
97
void ConjugateFMethod::update(RNG *rng)
98
    void ConjugateFMethod::update(RNG *rng)
98
{
99
    {
99
    /* There are three levels 
100
  /* There are three levels 
100
101
101
       X - variance parameter with F(m,1) prior
102
     X - variance parameter with F(m,1) prior
102
103
103
           X is represented internally using a redundant parameterization
104
           X is represented internally using a redundant parameterization
104
           X = tau/scale^2
105
           X = tau/scale^2
105
           where tau has a Chi-square distribution on m degrees of freedom
106
           where tau has a Chi-square distribution on m degrees of freedom
106
           and scale has a standard normal distribution
107
           and scale has a standard normal distribution
107
           
108
           
108
       Y - Stochastic children of X, which are normally distributed 
109
     Y - Stochastic children of X, which are normally distributed 
109
           random effects with precision that is a scale function of X
110
           random effects with precision that is a scale function of X
110
111
111
       Z - Stochastic Children of Y, which are normally distributed 
112
     Z - Stochastic Children of Y, which are normally distributed 
112
           random variables forming a linear model in X
113
           random variables forming a linear model in X
113
    */
114
  */
114
115
115
    /* Reparametrize from Y to scaled Y*/
116
  /* Reparametrize from Y to scaled Y*/
116
117
117
    vector<StochasticNode const*> const &Ynodes = _gv1->stochasticChildren();
118
  vector<StochasticNode const*> const &Ynodes = _gv1->stochasticChildren();
118
    unsigned int J = Ynodes.size();
119
  unsigned int J = Ynodes.size();
119
    vector<double> scaledY(J);
120
  vector<double> scaledY(J);
120
    for (unsigned int j = 0; j < J; ++j) {
121
  for (unsigned int j = 0; j < J; ++j) {
121
    scaledY[j] -= getValue(Ynodes[j], _chain) - getMean(Ynodes[j], _chain);
122
        scaledY[j] -= getValue(Ynodes[j], _chain) - getMean(Ynodes[j], _chain);
122
    scaledY[j] /= _scale; 
123
        scaledY[j] /= _scale; 
123
    }
124
  }
124
125
125
    /* 
126
  /* 
126
       Update _tau, which has a conjugate gamma distribution, conditional
127
     Update _tau, which has a conjugate gamma distribution, conditional
127
       on scaledY. (Note that _tau is independent of _scale). 
128
     on scaledY. (Note that _tau is independent of _scale). 
128
    */
129
  */
129
130
130
    //Prior
131
  //Prior
131
    vector<Node const *> const &Xparam = _gv1->nodes()[0]->parents();
132
  vector<Node const *> const &Xparam = _gv1->nodes()[0]->parents();
132
    double tau_shape = *Xparam[0]->value(_chain)/2; 
133
  double tau_shape = *Xparam[0]->value(_chain)/2; 
133
    double tau_rate = 0.5;
134
  double tau_rate = 0.5;
134
    
135
    
135
    //Likelihood 
136
  //Likelihood 
136
    double *coef = 0;
137
  double *coef = 0;
137
    bool empty = _gv1->deterministicChildren().empty();
138
  bool empty = _gv1->deterministicChildren().empty();
138
    if (!empty && _coef == 0) {
139
  if (!empty && _coef == 0) {
139
        coef = new double[J];
140
        coef = new double[J];
140
        calCoef(coef, _gv1, _chain);
141
        calCoef(coef, _gv1, _chain);
141
    }
142
  }
142
    else {
143
  else {
143
    coef = _coef;
144
        coef = _coef;
144
    }
145
  }
145
146
146
    for (unsigned int j = 0; j < J; ++j) {
147
  for (unsigned int j = 0; j < J; ++j) {
147
    double coef_j = empty ? 1 : coef[j];
148
        double coef_j = empty ? 1 : coef[j];
148
    if (coef_j > 0) {
149
        if (coef_j > 0) {
149
        tau_shape += 0.5;
150
      tau_shape += 0.5;
150
        tau_rate += coef_j * scaledY[j] * scaledY[j] / 2;
151
      tau_rate += coef_j * scaledY[j] * scaledY[j] / 2;
151
  }
152
    }
152
      }
153
  }
153
    if (coef != _coef) {
154
  if (coef != _coef) {
154
    delete [] coef;
155
        delete [] coef;
155
    coef = 0;
156
        coef = 0;
156
    }
157
  }
157
158
158
    // Sample from the posterior
159
  // Sample from the posterior
159
    _tau = rgamma(tau_shape, 1/tau_rate, rng);
160
  _tau = rgamma(tau_shape, 1/tau_rate, rng);
160
161
161
    /* 
162
  /* 
162
       Update _scale, which has a conjugate normal distribution given
163
     Update _scale, which has a conjugate normal distribution given
163
       Y and scaledY (Note that it is independent of _tau)
164
     Y and scaledY (Note that it is independent of _tau)
164
    */
165
  */
165
    vector<double> Yold(J); //Current value of Y
166
  vector<double> Yold(J); //Current value of Y
166
    vector<double> Ymean(J); //Prior mean of Y
167
  vector<double> Ymean(J); //Prior mean of Y
167
    for (unsigned int j = 0; j < J; ++j) {
168
  for (unsigned int j = 0; j < J; ++j) {
168
    Yold[j] = getValue(Ynodes[j], _chain);
169
        Yold[j] = getValue(Ynodes[j], _chain);
169
    Ymean[j] = getMean(Ynodes[j], _chain);
170
        Ymean[j] = getMean(Ynodes[j], _chain);
170
    }
171
  }
171
172
172
    vector<StochasticNode const *> const &Znodes = _gv2->stochasticChildren();
173
  vector<StochasticNode const *> const &Znodes = _gv2->stochasticChildren();
173
    unsigned int I = Znodes.size();
174
  unsigned int I = Znodes.size();
174
175
175
    /* 
176
  /* 
176
       Zprecision[i] is the precision of Z[i]
177
     Zprecision[i] is the precision of Z[i]
177
       Zresidual[i] is the difference between Z[i] and its prior mean
178
     Zresidual[i] is the difference between Z[i] and its prior mean
178
       when scaledY==0
179
     when scaledY==0
179
       Zlp[i] is the contribution of scaledY to the linear predictor
180
     Zlp[i] is the contribution of scaledY to the linear predictor
180
       (excluding intercept terms)
181
     (excluding intercept terms)
181
    */
182
  */
182
    vector<double> Zlp(I), Zresidual(I), Zprecision(I);
183
  vector<double> Zlp(I), Zresidual(I), Zprecision(I);
183
    for (unsigned int i = 0; i < I; ++i) {
184
  for (unsigned int i = 0; i < I; ++i) {
184
    Zresidual[i] = *Znodes[i]->value(_chain);
185
        Zresidual[i] = *Znodes[i]->value(_chain);
185
    Zlp[i] = getMean(Znodes[i], _chain);
186
        Zlp[i] = getMean(Znodes[i], _chain);
186
        Zprecision[i] = getScale(Znodes[i], _chain);
187
      Zprecision[i] = getScale(Znodes[i], _chain);
187
    }
188
  }
188
    _gv2->setValue(Ymean, _chain);
189
  _gv2->setValue(Ymean, _chain);
189
    for (unsigned int i = 0; i < I; ++i) {
190
  for (unsigned int i = 0; i < I; ++i) {
190
    Zresidual[i] -= getMean(Znodes[i], _chain);
191
        Zresidual[i] -= getMean(Znodes[i], _chain);
191
    Zlp[i] -= getMean(Znodes[i], _chain);
192
        Zlp[i] -= getMean(Znodes[i], _chain);
192
    Zlp[i] /= _scale; 
193
        Zlp[i] /= _scale; 
193
    }
194
  }
194
    _gv2->setValue(Yold, _chain);
195
  _gv2->setValue(Yold, _chain);
195
196
196
    double mean_scale = 0, precision_scale=1;
197
  double mean_scale = 0, precision_scale=1;
197
    for (unsigned int i = 0; i < I; ++i) {
198
  for (unsigned int i = 0; i < I; ++i) {
198
    precision_scale += Zlp[i] * Zlp[i] * Zprecision[i];
199
        precision_scale += Zlp[i] * Zlp[i] * Zprecision[i];
199
    mean_scale += Zresidual[i] * Zprecision[i] * Zlp[i];
200
        mean_scale += Zresidual[i] * Zprecision[i] * Zlp[i];
200
    }
201
  }
201
    mean_scale /= precision_scale;
202
  mean_scale /= precision_scale;
202
203
203
    for(unsigned int r = 0; r <= MAX_RNORM; r++) {
204
  for(unsigned int r = 0; r <= MAX_RNORM; r++) {
204
    if (r == MAX_RNORM) {
205
        if (r == MAX_RNORM) {
205
        throwRuntimeError("Degenerate scale in Conjugate F Sampler");
206
      throwRuntimeError("Degenerate scale in Conjugate F Sampler");
206
  }
207
      }
207
    _scale = rnorm(mean_scale, 1/sqrt(precision_scale), rng);
208
        _scale = rnorm(mean_scale, 1/sqrt(precision_scale), rng);
208
    if (fabs(_scale) > MIN_SCALE) {
209
        if (fabs(_scale) > MIN_SCALE) {
209
        break;
210
      break;
210
  }
211
    }
211
      }
212
  }
212
213
213
    double xnew = _tau/(_scale * _scale);
214
  double xnew = _tau/(_scale * _scale);
214
    //std::cout << _tau << " " << _scale << " " << xnew << std::endl;
215
  //std::cout << _tau << " " << _scale << " " << xnew << std::endl;
215
    _gv1->setValue(&xnew, 1, _chain);
216
  _gv1->setValue(&xnew, 1, _chain);
217
    }
218
216
}
219
}