[Rsiena-commits] r307 - in pkg: RSiena RSiena/R RSiena/src RSiena/src/model RSiena/src/model/effects RSiena/src/model/variables RSienaTest RSienaTest/doc RSienaTest/man RSienaTest/src RSienaTest/src/model RSienaTest/src/model/effects RSienaTest/src/model/variables

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Fri May 12 15:51:02 CEST 2017


Author: tomsnijders
Date: 2017-05-12 15:51:01 +0200 (Fri, 12 May 2017)
New Revision: 307

Added:
   pkg/RSiena/src/init.cpp
   pkg/RSienaTest/src/model/effects/BothDegreesEffect.cpp
   pkg/RSienaTest/src/model/effects/BothDegreesEffect.h
Modified:
   pkg/RSiena/ChangeLog
   pkg/RSiena/DESCRIPTION
   pkg/RSiena/R/initializeFRAN.r
   pkg/RSiena/src/model/StatisticCalculator.cpp
   pkg/RSiena/src/model/effects/EffectFactory.cpp
   pkg/RSiena/src/model/variables/BehaviorVariable.cpp
   pkg/RSiena/src/model/variables/BehaviorVariable.h
   pkg/RSienaTest/ChangeLog
   pkg/RSienaTest/DESCRIPTION
   pkg/RSienaTest/doc/Siena_algorithms.tex
   pkg/RSienaTest/man/RSiena-package.Rd
   pkg/RSienaTest/src/model/StatisticCalculator.cpp
   pkg/RSienaTest/src/model/effects/EffectFactory.cpp
   pkg/RSienaTest/src/model/variables/BehaviorVariable.cpp
   pkg/RSienaTest/src/model/variables/BehaviorVariable.h
   pkg/RSienaTest/src/siena07internals.cpp
Log:
Version 1.1-307, corrections mainly for behModelType=2

Modified: pkg/RSiena/ChangeLog
===================================================================
--- pkg/RSiena/ChangeLog	2017-05-10 15:32:04 UTC (rev 306)
+++ pkg/RSiena/ChangeLog	2017-05-12 13:51:01 UTC (rev 307)
@@ -4,6 +4,13 @@
      in the algorithm object. This required only a change in 
      NetworkVariable::checkAlterAgreement; documented in Siena_algorithms.
 
+2017-05-12 R-Forge Revision 307
+Changes in RSiena and RSienaTest:
+   * Operation of option 'absorb' (behModelType=2) corrected.
+Changes in RSienaTest:
+   * BothDegreesEffect.h and .cpp added to repository.
+
+
 2017-05-10 R-Forge Revision 306
 Changes in RSiena and RSienaTest:
    * Added updateSpecification (effectsMethods)

Modified: pkg/RSiena/DESCRIPTION
===================================================================
--- pkg/RSiena/DESCRIPTION	2017-05-10 15:32:04 UTC (rev 306)
+++ pkg/RSiena/DESCRIPTION	2017-05-12 13:51:01 UTC (rev 307)
@@ -1,8 +1,8 @@
 Package: RSiena
 Type: Package
 Title: Siena - Simulation Investigation for Empirical Network Analysis
-Version: 1.1-306
-Date: 2017-05-10
+Version: 1.1-307
+Date: 2017-05-12
 Author: Ruth Ripley, Krists Boitmanis, Tom A.B. Snijders, Felix Schoenenberger
 Depends: R (>= 2.15.0), utils
 Imports: Matrix, tcltk, lattice, parallel, MASS, methods

Modified: pkg/RSiena/R/initializeFRAN.r
===================================================================
--- pkg/RSiena/R/initializeFRAN.r	2017-05-10 15:32:04 UTC (rev 306)
+++ pkg/RSiena/R/initializeFRAN.r	2017-05-12 13:51:01 UTC (rev 307)
@@ -688,7 +688,6 @@
 		simpleRates <- FALSE
 	}
 	z$simpleRates <- simpleRates
-
 	ans <- .Call("setupModelOptions", PACKAGE=pkgname,
                  pData, pModel, MAXDEGREE, UNIVERSALOFFSET, CONDVAR, CONDTARGET,
                  profileData, z$parallelTesting, MODELTYPE, BEHMODELTYPE,

Added: pkg/RSiena/src/init.cpp
===================================================================
--- pkg/RSiena/src/init.cpp	                        (rev 0)
+++ pkg/RSiena/src/init.cpp	2017-05-12 13:51:01 UTC (rev 307)
@@ -0,0 +1,20 @@
+// Generated by the rstantools package
+
+
+#include <R.h>
+#include <Rinternals.h>
+#include <R_ext/Rdynload.h>
+#include <R_ext/Visibility.h>
+#include <Rversion.h>
+
+
+static const R_CallMethodDef CallEntries[] = {
+  {NULL, NULL, 0}
+};
+
+
+void attribute_visible R_init_RSiena(DllInfo *dll) {
+  // next line is necessary to avoid a NOTE from R CMD check
+  R_registerRoutines(dll, NULL, CallEntries, NULL, NULL);
+  R_useDynamicSymbols(dll, TRUE); // necessary for .onLoad() to work
+}


Property changes on: pkg/RSiena/src/init.cpp
___________________________________________________________________
Added: svn:eol-style
   + native

Modified: pkg/RSiena/src/model/StatisticCalculator.cpp
===================================================================
--- pkg/RSiena/src/model/StatisticCalculator.cpp	2017-05-10 15:32:04 UTC (rev 306)
+++ pkg/RSiena/src/model/StatisticCalculator.cpp	2017-05-12 13:51:01 UTC (rev 307)
@@ -652,8 +652,6 @@
 
 	double * currentValues  = new double[pBehaviorData->n()];
 
-	bool modelAbsorb = (BehaviorModelType(pBehaviorData->behModelType()) == ABSORB);
-
 	for (int i = 0; i < pBehaviorData->n(); i++)
 	{
 		currentValues[i] = currentState[i] - pBehaviorData->overallMean();
@@ -731,8 +729,8 @@
 				// no change gives no contribution
 				this->lstaticChangeContributions.at(pInfo).at(e)[1] = 0;
 				// calculate the contribution of downward change
-				if (((currentValues[e] > pBehaviorData->min()) &&
-						(!pBehaviorData->upOnly(this->lperiod))) || modelAbsorb)
+				if ((currentValues[e] > pBehaviorData->min()) &&
+						(!pBehaviorData->upOnly(this->lperiod)))
 				{
 					this->lstaticChangeContributions.at(pInfo).at(e)[0] =
 									pEffect->calculateChangeContribution(e,-1);
@@ -742,8 +740,8 @@
 					this->lstaticChangeContributions.at(pInfo).at(e)[0] = R_NaN;
 				}
 				// calculate the contribution of upward change
-				if (((currentValues[e] < pBehaviorData->max()) &&
-						(!pBehaviorData->downOnly(this->lperiod))) || modelAbsorb)
+				if ((currentValues[e] < pBehaviorData->max()) &&
+						(!pBehaviorData->downOnly(this->lperiod)))
 				{
 					this->lstaticChangeContributions.at(pInfo).at(e)[2] =
 									pEffect->calculateChangeContribution(e,1);
@@ -923,98 +921,101 @@
 
 	// for each covariate-defined setting, calculate
 	// setting*difference network and sum.
-	const vector<string> & rSettingNames = pNetworkData->rSettingNames();
+//	const vector<pair<string, string> > & rSettingNames =
+//			pNetworkData->rSettingNames();
 	//Rprintf("basic rate change %d size %d \n", pDifference->tieCount(),
 	//	rSettingNames.size());
-	for (unsigned i = 0; i < rSettingNames.size();
-		 i++)
-	{
-		if (!this->lsettingDistances[pNetworkData][rSettingNames[i]])
-		{
-			int * array =
-				new int[pNetworkData->observationCount() - 1];
+	// for (unsigned i = 0; i < rSettingNames.size();
+	// 	 i++)
+	// {
+	// 	if (!this->lsettingDistances[pNetworkData][rSettingNames[i]])
+	// 	{
+	// 		int * array =
+	// 			new int[pNetworkData->observationCount() - 1];
 
-			this->lsettingDistances[pNetworkData][rSettingNames[i]] = array;
-		}
-		// universal
-		int distance = pDifference->tieCount();
-		// primary
-		if (i == 1)
-		{
-			const Network * pNetwork =
-				this->lpState->pNetwork(pNetworkData->name());
-			// create a network representing primary settings
-			Network * settingNetwork = new
-				Network(pNetwork->n(), pNetwork->m());
-			for (int ego = 0; ego < pNetwork->n(); ego++)
-			{
-				vector<int> * setting = primarySetting(pNetwork, ego);
-				for (unsigned alter = 0; alter < setting->size(); alter++)
-				{
-					settingNetwork->setTieValue(ego, alter, 1);
-				}
-				delete(setting);
-			}
-			for (TieIterator iter = pDifference->ties();
-				 iter.valid();
-				 iter.next())
-			{
+	// 		this->lsettingDistances[pNetworkData][rSettingNames[i]] = array;
+	// 	}
+	// 	// universal
+	// 	int distance = pDifference->tieCount();
+	// 	// primary
+	// 	if (i == 1)
+	// 	{
+	// 		const Network * pNetwork =
+	// 			this->lpState->pNetwork(pNetworkData->name());
+	// 		// create a network representing primary settings
+	// 		Network * settingNetwork = new
+	// 			Network(pNetwork->n(), pNetwork->m());
+	// 		for (int ego = 0; ego < pNetwork->n(); ego++)
+	// 		{
+	// 			vector<int> * setting = primarySetting(pNetwork, ego);
+	// 			for (unsigned alter = 0; alter < setting->size(); alter++)
+	// 			{
+	// 				settingNetwork->setTieValue(ego, alter, 1);
+	// 			}
+	// 			delete(setting);
+	// 		}
+	// 		for (TieIterator iter = pDifference->ties();
+	// 			 iter.valid();
+	// 			 iter.next())
+	// 		{
 
-				{
-					if (settingNetwork->tieValue(iter.ego(),
-							iter.alter()) == 0)
-					{
-						distance--;
-					}
-				}
-			}
-			delete settingNetwork;
-		}
-		// for each covariate-defined setting, calculate
-		// setting*difference network and sum.
-		if (i > 1)
-		{
-			ConstantDyadicCovariate * pConstantDyadicCovariate =
-				this->lpData->pConstantDyadicCovariate(rSettingNames[i]);
-			ChangingDyadicCovariate * pChangingDyadicCovariate =
-				this->lpData->pChangingDyadicCovariate(rSettingNames[i]);
+	// 			{
+	// 				if (settingNetwork->tieValue(iter.ego(),
+	// 						iter.alter()) == 0)
+	// 				{
+	// 					distance--;
+	// 				}
+	// 			}
+	// 		}
+	// 		delete settingNetwork;
+	// 	}
+	// 	// for each covariate-defined setting, calculate
+	// 	// setting*difference network and sum.
+	// 	if (i > 1)
+	// 	{
+	// 		ConstantDyadicCovariate * pConstantDyadicCovariate =
+	// 			this->lpData->pConstantDyadicCovariate(rSettingNames[i]);
+	// 		ChangingDyadicCovariate * pChangingDyadicCovariate =
+	// 			this->lpData->pChangingDyadicCovariate(rSettingNames[i]);
 
-			for (TieIterator iter = pDifference->ties();
-				 iter.valid();
-				 iter.next())
-			{
+	// 		for (TieIterator iter = pDifference->ties();
+	// 			 iter.valid();
+	// 			 iter.next())
+	// 		{
 
-				if (pConstantDyadicCovariate)
-				{
-					if (pConstantDyadicCovariate->value(iter.ego(),
-							iter.alter()) == 0)
-					{
-						distance--;
-					}
-				}
-				else if (pChangingDyadicCovariate)
-				{
-					if (pChangingDyadicCovariate->value(iter.ego(),
-							iter.alter(),
-							this->lperiod) == 0)
-					{
-						distance--;
-					}
-				}
-				else
-				{
-					throw logic_error(
-						"No dyadic covariate named '" +
-						rSettingNames[i] +
-						"'.");
-				}
-			}
-		}
-		// store distance for later use
-		//	Rprintf(" cov %d %d\n", i, distance);
-		this->lsettingDistances[pNetworkData][rSettingNames[i]]
-			[this->lperiod] = distance;
-	}
+	// 			if (pConstantDyadicCovariate)
+	// 			{
+	// 				if (pConstantDyadicCovariate->value(iter.ego(),
+	// 						iter.alter()) == 0)
+	// 				{
+	// 					distance--;
+	// 				}
+	// 			}
+	// 			else if (pChangingDyadicCovariate)
+	// 			{
+	// 				if (pChangingDyadicCovariate->value(iter.ego(),
+	// 						iter.alter(),
+	// 						this->lperiod) == 0)
+	// 				{
+	// 					distance--;
+	// 				}
+	// 			}
+	// 			else
+	// 			{
+	// 				throw logic_error(
+	// 					"No dyadic covariate named '" +
+	// 					rSettingNames[i] +
+	// 					"'.");
+	// 			}
+	// 		}
+	// 	}
+	// 	// store distance for later use
+	// 	//	Rprintf(" cov %d %d\n", i, distance);
+	// 	this->lsettingDistances[pNetworkData][rSettingNames[i]]
+	// 		[this->lperiod] = distance;
+	// }
+	// the universal setting can explain everything
+//	calcDifferences(pNetworkData, pDifference);
 
 	// Loop through the rate effects, calculate the statistics,
 	// and store them.

Modified: pkg/RSiena/src/model/effects/EffectFactory.cpp
===================================================================
--- pkg/RSiena/src/model/effects/EffectFactory.cpp	2017-05-10 15:32:04 UTC (rev 306)
+++ pkg/RSiena/src/model/effects/EffectFactory.cpp	2017-05-12 13:51:01 UTC (rev 307)
@@ -525,7 +525,7 @@
 	{
 		pEffect = new FourCyclesEffect(pEffectInfo, false);
 	}
-	else if ((effectName == "gwespFF")|(effectName == "gwesp"))
+	else if ((effectName == "gwespFF")||(effectName == "gwesp"))
 	{
 		EgocentricConfigurationTable * (NetworkCache::*mytable)() const =
 			&NetworkCache::pTwoPathTable;
@@ -715,7 +715,7 @@
  		pEffect = new GenericNetworkEffect(pEffectInfo,
 			pFunction);
 	}
-	else if ((effectName == "gwespFBMix")|(effectName == "gwespMix"))
+	else if ((effectName == "gwespFBMix")||(effectName == "gwespMix"))
  	{
 		EgocentricConfigurationTable * (NetworkCache::*mytable)() const =
 			&NetworkCache::pInStarTable;

Modified: pkg/RSiena/src/model/variables/BehaviorVariable.cpp
===================================================================
--- pkg/RSiena/src/model/variables/BehaviorVariable.cpp	2017-05-10 15:32:04 UTC (rev 306)
+++ pkg/RSiena/src/model/variables/BehaviorVariable.cpp	2017-05-12 13:51:01 UTC (rev 307)
@@ -50,7 +50,6 @@
 	this->lendowmentEffectContribution = new double * [3];
 	this->lcreationEffectContribution = new double * [3];
 	this->lprobabilities = new double[3];
-	this->lqrobabilities = new double[3];
 
 	for (int i = 0; i < 3; i++)
 	{
@@ -61,11 +60,9 @@
 		this->lcreationEffectContribution[i] =
 			new double[pSimulation->pModel()->rCreationEffects(pData->name()).size()];
 		this->lprobabilities[i] = 0;
-		this->lqrobabilities[i] = 0;
 	}
 
 	this->lbehaviorModelType = BehaviorModelType(pData->behModelType());
-	this->lmodelAbsorb = (this->lbehaviorModelType == ABSORB);
 }
 
 
@@ -79,7 +76,6 @@
 	this->lpData = 0;
 	this->lvalues = 0;
 	delete[] this->lprobabilities;
-	delete[] this->lqrobabilities;
 	// Delete arrays of contributions
 
 	for (int i = 0; i < 3; i++)
@@ -97,7 +93,6 @@
 	this->lendowmentEffectContribution = 0;
 	this->lcreationEffectContribution = 0;
 	this->lprobabilities = 0;
-	this->lqrobabilities = 0;
 
 	// no need to delete lpChangeContribution since this is
 	// handled by the MiniStep
@@ -225,6 +220,17 @@
 	return this->lbehaviorModelType;
 }
 
+/**
+ * Returns whether the model type is absorb.
+ */
+
+bool BehaviorVariable::behaviorModelTypeABSORB() const
+{
+	return (this->lbehaviorModelType == ABSORB);
+}
+
+
+
 // ----------------------------------------------------------------------------
 // Section: Initialization at the beginning of a period
 // ----------------------------------------------------------------------------
@@ -242,6 +248,8 @@
 	{
 		this->lvalues[i] = this->lpData->value(period, i);
 	}
+
+	this->behaviorModelType(this->lpData->behModelType());
 }
 
 
@@ -285,7 +293,7 @@
 
 	if (this->pSimulation()->pModel()->needScores())
 	{
-		this->accumulateScores(this->lvalues[actor], difference + 1,
+		this->accumulateScores(difference + 1,
 			this->lupPossible,
 			this->ldownPossible);
 	}
@@ -341,6 +349,8 @@
 	this->lupPossible = true;
 	this->ldownPossible = true;
 	int currentValue = this->lvalues[actor];
+	bool ismax = (currentValue >= this->lpData->max());
+	bool ismin = (currentValue <= this->lpData->min());
 
 	int evaluationEffectCount = this->pEvaluationFunction()->rEffects().size();
 	int endowmentEffectCount = this->pEndowmentFunction()->rEffects().size();
@@ -386,19 +396,10 @@
 		}
 	}
 
-	// Calculate the probability for downward change.
+	// Calculate the objective function for downward change.
 	// Defer exp until we can subtract the largest to avoid overflow.
-
-	if (((currentValue > this->lpData->min()) &&
-		(!this->lpData->upOnly(this->period()))) || (this->lmodelAbsorb))
+	if (ismin || this->lpData->upOnly(this->period()))
 	{
-		this->lprobabilities[0] =
-			this->totalEvaluationContribution(actor, -1) +
-				this->totalEndowmentContribution(actor, -1);
-		maxValue = max(maxValue, this->lprobabilities[0]);
-	}
-	else
-	{
 		this->lprobabilities[0] = 0;
 		this->ldownPossible = false;
 		for (unsigned i = 0; i < pEvaluationFunction()->rEffects().size(); i++)
@@ -417,22 +418,20 @@
 			this->lcreationEffectContribution[0][i] = R_NaN;
 		}
 	}
+	else
+	{
+		this->lprobabilities[0] =
+			this->totalEvaluationContribution(actor, -1) +
+				this->totalEndowmentContribution(actor, -1);
+		maxValue = max(maxValue, this->lprobabilities[0]);
+	}
 
 	// No change means zero contribution, but exp(0) = 1
 	this->lprobabilities[1] = 0;
 
-	// Calculate the probability for upward change
-
-	if (((currentValue < this->lpData->max()) &&
-		(!this->lpData->downOnly(this->period())))  || (this->lmodelAbsorb))
+	// Calculate the objective function for upward change
+	if (ismax || this->lpData->downOnly(this->period()))
 	{
-		this->lprobabilities[2] =
-			this->totalEvaluationContribution(actor, 1) +
-				this->totalCreationContribution(actor, 1);
-		maxValue = max(maxValue, this->lprobabilities[2]);
-	}
-	else
-	{
 		this->lprobabilities[2] = 0;
 		this->lupPossible = false;
 			for (unsigned i = 0; i < pEvaluationFunction()->rEffects().size(); i++)
@@ -450,42 +449,48 @@
 			this->lcreationEffectContribution[2][i] = R_NaN;
 		}
 	}
+	else
+	{
+		this->lprobabilities[2] =
+			this->totalEvaluationContribution(actor, 1) +
+				this->totalCreationContribution(actor, 1);
+		maxValue = max(maxValue, this->lprobabilities[2]);
+	}
 
+		// turn objective function values into probabilities
 	if (this->ldownPossible)
 	{
-		this->lprobabilities[0] -= maxValue;
-		this->lprobabilities[0] = exp(this->lprobabilities[0]);
+		this->lprobabilities[0] = exp(this->lprobabilities[0] - maxValue);
 	}
 	if (this->lupPossible)
 	{
-		this->lprobabilities[2] -= maxValue;
-		this->lprobabilities[2] = exp(this->lprobabilities[2]);
+		this->lprobabilities[2] = exp(this->lprobabilities[2] - maxValue);
 	}
 	this->lprobabilities[1]  = exp(-maxValue);
 
-	double sum = this->lprobabilities[0] + this->lprobabilities[1] +
-		this->lprobabilities[2];
-	this->lprobabilities[0] /= sum;
-	this->lprobabilities[1] /= sum;
-	this->lprobabilities[2] /= sum;
-
-	if (this->lmodelAbsorb)
+	double sum = 0;
+	if ((this->behaviorModelTypeABSORB()) && (ismin || ismax))
 	{
-		this->lqrobabilities = this->lprobabilities;
-		if (currentValue <= this->lpData->min())
+		if (ismin)
 		{
-			this->lprobabilities[1] =
-				this->lprobabilities[0] + this->lprobabilities[1];
-			this->lprobabilities[0] = 0;
-			this->ldownPossible = false;
+			sum = 2*this->lprobabilities[1] + this->lprobabilities[2];
+			this->lprobabilities[1] = 2*this->lprobabilities[1]/sum;
+			this->lprobabilities[2] /= sum;
 		}
-		if (currentValue >= this->lpData->max())
+		else
 		{
-			this->lprobabilities[1] =
-				this->lprobabilities[1] + this->lprobabilities[2];
-			this->lprobabilities[2] = 0;
-			this->lupPossible = false;
+			sum = 2*this->lprobabilities[1] + this->lprobabilities[0];
+			this->lprobabilities[1] = 2*this->lprobabilities[1]/sum;
+			this->lprobabilities[0] /= sum;
 		}
+		}
+	else
+	{
+		sum = this->lprobabilities[0] + this->lprobabilities[1] +
+							this->lprobabilities[2];
+		this->lprobabilities[0] /= sum;
+		this->lprobabilities[1] /= sum;
+		this->lprobabilities[2] /= sum;
 	}
 }
 
@@ -580,43 +585,27 @@
  * This function is called with arguments difference = 0, 1, 2.
  * 1 means no change.
  */
-void BehaviorVariable::accumulateScores(int currentVal, int difference,
+void BehaviorVariable::accumulateScores(int difference,
 	bool upPossible, bool downPossible) const
 {
 	for (unsigned i = 0;
 		i < this->pEvaluationFunction()->rEffects().size();
 		i++)
 	{
+// 		if (difference == 1) no change, but not initialised
+// 		{
+// 			this->levaluationEffectContribution[difference][i] = 0;
+// 		}
 		Effect * pEffect = this->pEvaluationFunction()->rEffects()[i];
-
-		double score = 0;
-
-		if ((difference == 1) && (this->lmodelAbsorb))
+		double score = this->levaluationEffectContribution[difference][i];
+		if (upPossible)
 		{
-			if (currentVal <= this->lpData->min())
-			{
-				score = this->levaluationEffectContribution[0][i] *
-							(this->lqrobabilities[0]/this->lprobabilities[1]);
-			}
-			if (currentVal >= this->lpData->max())
-			{
-				score = this->levaluationEffectContribution[2][i] *
-							(this->lqrobabilities[2]/this->lprobabilities[1]);
-			}
-		}
-		else if (difference != 1)
-		{
-			score = this->levaluationEffectContribution[difference][i];
-// which is 0 if difference==1, i.e., no change
-		}
-		if ((upPossible) || (this->lmodelAbsorb))
-		{
 			score -=
 				this->levaluationEffectContribution[2][i] *
 				this->lprobabilities[2];
 		}
 
-		if ((downPossible) || (this->lmodelAbsorb))
+		if (downPossible)
 		{
 			score -=
 				this->levaluationEffectContribution[0][i] *
@@ -652,7 +641,7 @@
 			score = this->lendowmentEffectContribution[difference][i];
 		}
 
-		if ((downPossible) || (this->lmodelAbsorb))
+		if (downPossible)
 		{
 			score -=
 				this->lendowmentEffectContribution[0][i] *
@@ -681,7 +670,7 @@
 			score = this->lcreationEffectContribution[difference][i];
 		}
 
-		if ((upPossible) || (this->lmodelAbsorb))
+		if (upPossible)
 		{
 			score -=
 				this->lcreationEffectContribution[2][i] *
@@ -752,8 +741,7 @@
 	this->calculateProbabilities(pMiniStep->ego());
 	if (this->pSimulation()->pModel()->needScores())
 	{
-		this->accumulateScores(this->lvalues[pMiniStep->ego()],
-			pBehaviorChange->difference() + 1,
+		this->accumulateScores(pBehaviorChange->difference() + 1,
 				this->lupPossible, this->ldownPossible);
 	}
 	if (this->pSimulation()->pModel()->needDerivatives())

Modified: pkg/RSiena/src/model/variables/BehaviorVariable.h
===================================================================
--- pkg/RSiena/src/model/variables/BehaviorVariable.h	2017-05-10 15:32:04 UTC (rev 306)
+++ pkg/RSiena/src/model/variables/BehaviorVariable.h	2017-05-12 13:51:01 UTC (rev 307)
@@ -49,6 +49,7 @@
 
 	void behaviorModelType(int type);
 	virtual BehaviorModelType behaviorModelType() const;
+	virtual bool behaviorModelTypeABSORB() const;
 
 	virtual void makeChange(int actor);
 
@@ -78,7 +79,7 @@
 		int difference) const;
 	double totalCreationContribution(int actor,
 		int difference) const;
-	void accumulateScores(int currentVal, int difference, bool UpPossible,
+	void accumulateScores(int difference, bool UpPossible,
 		bool downPossible) const;
 	void calculateProbabilities(int actor);
 	void accumulateDerivatives() const;
@@ -114,9 +115,6 @@
 
 	double * lprobabilities;
 
-	// Selection probability per each difference, for the behModelType ABSORB
-	double * lqrobabilities;
-
 	// Indicates if upward change is possible in the current situation
 	bool lupPossible;
 
@@ -128,9 +126,6 @@
 
 	// the model type
 	BehaviorModelType lbehaviorModelType;
-	
-	// whether the model type is ABSORB
-	bool lmodelAbsorb;
 };
 
 }

Modified: pkg/RSienaTest/ChangeLog
===================================================================
--- pkg/RSienaTest/ChangeLog	2017-05-10 15:32:04 UTC (rev 306)
+++ pkg/RSienaTest/ChangeLog	2017-05-12 13:51:01 UTC (rev 307)
@@ -4,6 +4,13 @@
      in the algorithm object. This required only a change in 
      NetworkVariable::checkAlterAgreement; documented in Siena_algorithms.
 
+2017-05-12 R-Forge Revision 307
+Changes in RSiena and RSienaTest:
+   * Operation of option 'absorb' (behModelType=2) corrected.
+Changes in RSienaTest:
+   * BothDegreesEffect.h and .cpp added to repository.
+
+
 2017-05-10 R-Forge Revision 306
 Changes in RSiena and RSienaTest:
    * Added updateSpecification (effectsMethods)

Modified: pkg/RSienaTest/DESCRIPTION
===================================================================
--- pkg/RSienaTest/DESCRIPTION	2017-05-10 15:32:04 UTC (rev 306)
+++ pkg/RSienaTest/DESCRIPTION	2017-05-12 13:51:01 UTC (rev 307)
@@ -1,8 +1,8 @@
 Package: RSienaTest
 Type: Package
 Title: Siena - Simulation Investigation for Empirical Network Analysis
-Version: 1.1-306
-Date: 2017-05-10
+Version: 1.1-307
+Date: 2017-05-12
 Author: Ruth Ripley, Krists Boitmanis, Tom A.B. Snijders, Felix Schoenenberger
 Depends: R (>= 2.15.0), utils
 Imports: Matrix, tcltk, lattice, parallel, MASS, RUnit, methods

Modified: pkg/RSienaTest/doc/Siena_algorithms.tex
===================================================================
--- pkg/RSienaTest/doc/Siena_algorithms.tex	2017-05-10 15:32:04 UTC (rev 306)
+++ pkg/RSienaTest/doc/Siena_algorithms.tex	2017-05-12 13:51:01 UTC (rev 307)
@@ -5428,7 +5428,7 @@
 The network model type then is used further in
 \texttt{NetworkVariable.cpp}.
 In  \texttt{BehaviorVariable.cpp}, the behavioral model type
-then is used further in its local form \texttt{lmodelAbsorb}.
+then is used further through the function \texttt{behaviorModelTypeABSORB}.
 
 \subsection{Behavior micro-step}
 Whenever actor $i$ may make a change in variable $Z$,
@@ -5482,8 +5482,8 @@
 Similarly, if $z_i = z^+ = \max\{\text{range}(Z)\}$, the probabilities are
 
 \[
-p_{i}(z^+ ; \beta, z^+, x) = \frac{\exp\big(f(i,z^+ )\big)  )\big)}
-             {\exp\big(f(i,z^+ - 1)\big) + \exp\big(f(i,z^+ )\big) \big)}
+p_{i}(z^+ ; \beta, z^+, x) = \frac{\exp\big(f(i,z^+ )\big)  }
+             {\exp\big(f(i,z^+ - 1)\big) + \exp\big(f(i,z^+ )\big) }
 \]
 and
 \[
@@ -5495,29 +5495,27 @@
 if $z_i = z^- = \min\{\text{range}(Z)\}$,
 the probabilities are
 \[
-p_{i}(z^- ; \beta, z^-, x) = \frac{\exp\big(f(i,z^- - 1)\big) + \exp\big(f(i,z^- )\big)}
-             {\exp\big(f(i,z^- - 1)\big) +
-             \exp\big(f(i,z^- )\big) + \exp\big(f(i,z^- + 1)\big)}
+p_{i}(z^- ; \beta, z^-, x) = \frac{2\, \exp\big(f(i,z^- )\big)}
+             {2\, \exp\big(f(i,z^- )\big) + \exp\big(f(i,z^- + 1)\big)}
 \]
 and
 \[
 p_{}(z^- + 1; \beta, z^-, x) = \frac{\exp\big(f(i,z^- + 1)\big) }
-             {\exp\big(f(i,z^- - 1)\big) +
-             \exp\big(f(i,z^- )\big) + \exp\big(f(i,z^- + 1)\big)} \ .
+             {2\, \exp\big(f(i,z^- )\big) + \exp\big(f(i,z^- + 1)\big)} \ .
 \]
 
 Similarly, if $z_i = z^+ = \max\{\text{range}(Z)\}$, the probabilities are
 
 \[
-p_{i}(z^+ ; \beta, z^+, x) = \frac{\exp\big(f(i,z^+ )\big) + \exp\big(f(i,z^+ +1 )\big)}
+p_{i}(z^+ ; \beta, z^+, x) = \frac{2\,\exp\big(f(i,z^+ )\big) }
              {\exp\big(f(i,z^+ - 1)\big)
-             + \exp\big(f(i,z^+ )\big) + \exp\big(f(i,z^+ + 1)\big)}
+             + 2\,\exp\big(f(i,z^+ )\big) }
 \]
 and
 \[
 p_{}(z^+ - 1; \beta, z^+, x) = \frac{\exp\big(f(i,z^+ - 1)\big) }
              {\exp\big(f(i,z^+ - 1)\big)
-             + \exp\big(f(i,z^+ )\big) + \exp\big(f(i,z^+ + 1)\big)} \ .
+             + 2\,\exp\big(f(i,z^+ )\big) } \ .
 \]
 
 \subsection{Score function}
@@ -5555,14 +5553,11 @@
 \end{equation}
 where
 \begin{equation}
-    \overline{\Delta_{ik}(. , z)} \,=\, \sum_{d=-1}^1  q_i(d) \,  \Delta_{ik}(d, z) \
+    \overline{\Delta_{ik}(. , z)} \,=\, \sum_{d=-1}^1  p_i(d) \,  \Delta_{ik}(d, z) \ .
 \end{equation}
 \end{subequations}
-for $q_i(d) = p_i(d)$.
 
 
-
-
 In the standard model (`\texttt{RESTRICT}'), for the boundary cases:\\
 if the current state is at the minimum, we have
 \begin{subequations}\label{p.r}
@@ -5588,10 +5583,9 @@
 \end{subequations}
 
 The scores in the boundary cases still are given by (\ref{sc1})
-but with $q_i(d) = \pi_i(d)$ given in (\ref{p.r}).
+but with $p_i(d) = \pi_i(d)$ given in (\ref{p.r}).
 
 
-
 In the new model (`\texttt{ABSORB}'), the probabilities are:
 
 \noindent
@@ -5614,67 +5608,26 @@
 \end{equation}
 \end{subequations}
 
+To calculate the scores in the second model type (`\texttt{ABSORB}'),
+for the boundary cases,
+we may note that this is based on a multinomial regression model
+with three options $\{-1, 0, 1\}$, of which the first two outcomes
+are combined. 
+Consider the case for the right boundary.
+The first two outcomes have the same value 
+$\Delta_{ik}(-1,z) = \Delta_{ik}(0,z) = 0$; the value is 0
+because this option means no change.
+A score function is a function of the sufficient statistic,
+and for this 3-option statistical model the sufficient statistic corresponds
+to the partition of the outcome space into  $\big\{\{-1, 0\}, \{1\} \big\}$.
+Therefore the score function that we need is again (\ref{sc1}),
+for the original probabilities $p_i$. However, since 
+$\Delta_{ik}(-1,z) = \Delta_{ik}(0,z) = 0$ and $\pi_i(1) = p_i(1)$,
+it does not matter whether we calculate (\ref{sc1})
+for $p_i$ or for $\pi_i$.
+ 
 
-In the second model type (`\texttt{ABSORB}'), however, the scores are different
-because we no longer are employing a multinomial regression model.
 
-We use the rule that, for a joint distribution of $(X,V)$
-with joint score function $J_{XV}$ and marginal score function $J_X$,
-\begin{equation}\label{s.xv}
-  J_X \,=\, \text{E}\big\{ J_{XV} \mid X \big\} \ .
-\end{equation}
-In  this case $V$ is the variable that takes all three values --1, 0, and +1,
- and $X$ takes the two values absorbed to the range of $Z$.
-
-So we have to use the conditional expectation of the score function
-for the three-outcome case.
-
-
-Therefore, if the current state is at the minimum,
-the score is
-\begin{subequations}\label{sc.a}
-\begin{equation}
-  J_k(0) \,=\, \frac{p_i(-1)}{\pi_i(0)}
-                  \Delta_{ik}(-1, z) \, - \, \overline{ \Delta_{ik}(. , z)} \ ,
-\hspace{1em}
-   J_k(1) \,=\,  \Delta_{ik}(1, z) \, - \, \overline{ \Delta_{ik}(. , z)}
-\end{equation}
-with $ \overline{ \Delta_{ik}(. , z)}$ as in (\ref{sc1}b),
-for $q_i(d) = p_i(d)$.
-If the current state is at the maximum,
-the score is
-\begin{equation}
-\hspace{1em}
-   J_k(-1) \,=\,  \Delta_{ik}(-1, z) \, - \, \overline{ \Delta_{ik}(. , z)}   \ ,
-\hspace{1em}
-  J_k(0) \,=\, \frac{p_i(1)}{\pi_i(0)}
-                  \Delta_{ik}(1, z) \, - \, \overline{ \Delta_{ik}(. , z)}
-\end{equation}
-\end{subequations}
-again with $ \overline{ \Delta_{ik}(. , z)}$ as in (\ref{sc1}b),
-for $q_i(d) = p_i(d)$.
-
-
-It can be checked that the score function for the outcome that is
-observed like it is originally remains identical;
-and that the expected value of the score function is 0.
-
-For arranging the computations, it may be noted that
-we still can use formulas (\ref{sc1}) but with $ \Delta_{ik}(0, z)$
-replaced by
-\begin{subequations}\label{sca2}
-\begin{equation}
-  \tilde \Delta(0) \,=\, \frac{p_i(-1) \, \Delta_{ik}(-1, z) } {p_i(-1) \,+\, p_i(0)} \
-\end{equation}
-for the lower boundary case, and
-\begin{equation}
-  \tilde \Delta(0) \,=\, \frac{p_i(1)\,  \Delta_{ik}(1, z) } {p_i(0) \,+\, p_i(1)} \
-\end{equation}
-\end{subequations}
-for the upper boundary case; and $q_i(d) = p_i(d)$ in (\ref{sc1}b).
-
-
-
 \newpage
 \begin{appendices}
 

Modified: pkg/RSienaTest/man/RSiena-package.Rd
===================================================================
--- pkg/RSienaTest/man/RSiena-package.Rd	2017-05-10 15:32:04 UTC (rev 306)
+++ pkg/RSienaTest/man/RSiena-package.Rd	2017-05-12 13:51:01 UTC (rev 307)
@@ -46,8 +46,8 @@
 \tabular{ll}{
 Package: \tab RSienaTest\cr
 Type: \tab Package\cr
-Version: \tab 1.1-306\cr
-Date: \tab 2017-05-10\cr
+Version: \tab 1.1-307\cr
+Date: \tab 2017-05-12\cr
 Depends: \tab R (>= 3.0.0)\cr
 Imports: \tab Matrix\cr
 Suggests: \tab tcltk, network, codetools, lattice, MASS, parallel,

Modified: pkg/RSienaTest/src/model/StatisticCalculator.cpp
===================================================================
--- pkg/RSienaTest/src/model/StatisticCalculator.cpp	2017-05-10 15:32:04 UTC (rev 306)
+++ pkg/RSienaTest/src/model/StatisticCalculator.cpp	2017-05-12 13:51:01 UTC (rev 307)
@@ -97,7 +97,7 @@
 	this->lpStateLessMissingsEtc = new State();
 	this->lneedActorStatistics = returnActorStatistics;
 	this->lcountStaticChangeContributions = 0;
-	
+
 	this->calculateStatistics();
 }
 
@@ -120,7 +120,7 @@
 	this->lpStateLessMissingsEtc = new State();
 	this->lneedActorStatistics = returnActorStatistics;
 	this->lcountStaticChangeContributions = returnStaticChangeContributions;
-	
+
 	this->calculateStatistics();
 }
 
@@ -426,8 +426,6 @@
 			pBehaviorData->name());
 	double* currentValues = new double[pBehaviorData->n()];
 
-	bool modelAbsorb = (BehaviorModelType(pBehaviorData->behModelType()) == ABSORB);
-
 	for (int i = 0; i < pBehaviorData->n(); i++) {
 		currentValues[i] = currentState[i] - pBehaviorData->overallMean();
 		if (pBehaviorData->missing(this->lperiod, i)
@@ -476,8 +474,8 @@
 				// no change gives no contribution
 				this->lstaticChangeContributions.at(pInfo).at(e)[1] = 0;
 				// calculate the contribution of downward change
-				if (((currentValues[e] > pBehaviorData->min())
-						&& (!pBehaviorData->upOnly(this->lperiod))) || modelAbsorb)
+				if ((currentValues[e] > pBehaviorData->min())
+						&& (!pBehaviorData->upOnly(this->lperiod)))
 				{
 					this->lstaticChangeContributions.at(pInfo).at(e)[0] =
 							pEffect->calculateChangeContribution(e, -1);
@@ -485,8 +483,8 @@
 					this->lstaticChangeContributions.at(pInfo).at(e)[0] = R_NaN;
 				}
 				// calculate the contribution of upward change
-				if (((currentValues[e] < pBehaviorData->max())
-						&& (!pBehaviorData->downOnly(this->lperiod))) || modelAbsorb)
+				if ((currentValues[e] < pBehaviorData->max())
+						&& (!pBehaviorData->downOnly(this->lperiod)))
 					{
 					this->lstaticChangeContributions.at(pInfo).at(e)[2] =
 							pEffect->calculateChangeContribution(e, 1);
@@ -820,7 +818,7 @@
 //					}
 //				}
 //			}
-			
+
 			delete pEffect;
 		}
 
@@ -845,8 +843,6 @@
 
 	double * currentValues  = new double[pBehaviorData->n()];
 
-	bool modelAbsorb = (BehaviorModelType(pBehaviorData->behModelType()) == ABSORB);
-
 	for (int i = 0; i < pBehaviorData->n(); i++)
 	{
 		currentValues[i] = currentState[i] - pBehaviorData->overallMean();
@@ -923,8 +919,8 @@
 				// no change gives no contribution
 				this->lstaticChangeContributions.at(pInfo).at(e)[1] = 0;
 				// calculate the contribution of downward change
-				if (((currentValues[e] > pBehaviorData->min()) &&
-						(!pBehaviorData->upOnly(this->lperiod))) || modelAbsorb)
+				if ((currentValues[e] > pBehaviorData->min()) &&
+						(!pBehaviorData->upOnly(this->lperiod)))
 				{
 					this->lstaticChangeContributions.at(pInfo).at(e)[0] =
 									pEffect->calculateChangeContribution(e,-1);
@@ -934,8 +930,8 @@
 					this->lstaticChangeContributions.at(pInfo).at(e)[0] = R_NaN;
 				}
 				// calculate the contribution of upward change
-				if (((currentValues[e] < pBehaviorData->max()) &&
-						(!pBehaviorData->downOnly(this->lperiod))) || modelAbsorb)
+				if ((currentValues[e] < pBehaviorData->max()) &&
+						(!pBehaviorData->downOnly(this->lperiod)))
 				{
 					this->lstaticChangeContributions.at(pInfo).at(e)[2] =
 									pEffect->calculateChangeContribution(e,1);
@@ -946,7 +942,7 @@
 				}
 			}
 		}
-		
+
 		delete pEffect;
 	}
 

Added: pkg/RSienaTest/src/model/effects/BothDegreesEffect.cpp
===================================================================
--- pkg/RSienaTest/src/model/effects/BothDegreesEffect.cpp	                        (rev 0)
+++ pkg/RSienaTest/src/model/effects/BothDegreesEffect.cpp	2017-05-12 13:51:01 UTC (rev 307)
@@ -0,0 +1,110 @@
+/******************************************************************************
+ * SIENA: Simulation Investigation for Empirical Network Analysis
+ *
+ * Web: http://www.stats.ox.ac.uk/~snijders/siena/
+ *
+ * File: BothDegreesEffect.cpp
+ *
+ * Description: This file contains the implementation of the
+ * BothDegreesEffect class.
+ *****************************************************************************/
+
+#include "BothDegreesEffect.h"
+#include "utils/SqrtTable.h"
+#include "network/Network.h"
+#include "model/variables/NetworkVariable.h"
+#include "model/EffectInfo.h"
+#include "Effect.h"
+
+namespace siena
+{
+
+/**
+ * Constructor.
+ */
+BothDegreesEffect::BothDegreesEffect(
+	const EffectInfo * pEffectInfo) : NetworkEffect(pEffectInfo)
+{
+	this->lroot = (pEffectInfo->internalEffectParameter() >= 2);
+	this->lsqrtTable = SqrtTable::instance();
+}
+
+/**
+ * Calculates the contribution of a tie flip to the given actor.
+ */
[TRUNCATED]

To get the complete diff run:
    svnlook diff /svnroot/rsiena -r 307


More information about the Rsiena-commits mailing list