marca doc NCD

This example looks at transient probabilites calculated in section 7 of http://www4.ncsu.edu/~billy/MARCA/marca_doc.ps.

input

  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
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
"""
The NCD example from marca_doc.ps.

"""
from __future__ import print_function, division, absolute_import

import numpy as np

import jsonctmctree.interface


def gen_discrete_simplex(n, k):
    """
    n : the population size
    k : the number of variables
    yield suffixes
    """
    if not n:
        yield [0] * k
    elif k == 1:
        yield [n]
    else:
        for i in range(n+1):
            prefix = [i]
            for suffix in gen_discrete_simplex(n - i, k - 1):
                yield prefix + suffix


def gen_triples(n):
    # initial state, final state, rate
    states = list(gen_discrete_simplex(n, 4))
    for sa in states:
        CPU, SM, FD, TERM = sa
        if CPU:
            SM_rate = 100 * ((CPU + SM + FD) / 128)**1.5
            if SM < 10:
                yield sa, [CPU-1, SM+1, FD, TERM], SM_rate
            if FD < 10:
                yield sa, [CPU-1, SM, FD+1, TERM], 0.05
            if TERM < 10:
                yield sa, [CPU-1, SM, FD, TERM+1], 0.002
        if CPU < 10:
            if SM:
                yield sa, [CPU+1, SM-1, FD, TERM], 0.2
            if FD:
                yield sa, [CPU+1, SM, FD-1, TERM], 1/30
            if TERM:
                yield sa, [CPU+1, SM, FD, TERM-1], 0.0001 * TERM


def get_process_definition():
    triples = list(gen_triples(10))
    row_states, column_states, rates = zip(*triples)
    return dict(
            row_states = list(row_states),
            column_states = list(column_states),
            transition_rates = list(rates))


def get_scene(t, observations):
    # t is the time.

    tree = dict(
            row_nodes = [0],
            column_nodes = [1],
            edge_rate_scaling_factors = [t],
            edge_processes = [0])

    root_prior = dict(
            states = [
                [4, 0, 0, 6],
                [4, 1, 0, 5],
                [4, 0, 1, 5],
                [4, 2, 0, 4]],
            probabilities = [
                0.25, 0.25, 0.25, 0.25])

    # The first three observations are specifically requested.
    # The sum of the last four observations gives the total
    observed_data = dict(
            nodes = [1, 1, 1, 1],
            variables = [0, 1, 2, 3],
            iid_observations = observations)

    process_definition = get_process_definition()

    # Assemble the scene.
    scene = dict(
            node_count = 2,
            process_count = 1,
            state_space_shape = [11, 11, 11, 11],
            tree = tree,
            root_prior = root_prior,
            process_definitions = [process_definition],
            observed_data = observed_data)
    
    return scene


def main():

    # Define some requests.
    # These include the log likelihood
    # and the sum of transition count expectations.
    log_likelihood_request = dict(property = 'DNNLOGL')

    # Just get likelihoods for all possible observations.
    observations = list(gen_discrete_simplex(10, 4))
    obs_to_idx = dict((tuple(x), i) for i, x in enumerate(observations))

    for t in 10, 25:
        print('t:', t)
        print('-----')
        print()

        # Get the per-state log likelihoods.
        scene = get_scene(t, observations)
        j_in = {
                "scene" : scene,
                "requests" : [log_likelihood_request],
                }
        j_out = jsonctmctree.interface.process_json_in(j_in)
        log_likelihoods = j_out['responses'][0]
        likelihoods = np.exp(log_likelihoods)

        print('custom state probability requests:')
        for obs in (
                (0, 0, 0, 10),
                (6, 0, 0, 4),
                (4, 0, 0, 6)):
            p = likelihoods[obs_to_idx[obs]]
            print(list(obs), p, sep='\t')
        print()

        print('cumulative probabilities for constrained states:')
        cumulative = 0
        for obs in (
                (4, 0, 0, 6),
                (4, 1, 0, 5),
                (4, 0, 1, 5),
                (4, 2, 0, 4),
                (4, 1, 1, 4),
                (4, 0, 2, 4),
                (4, 3, 0, 3),
                (4, 2, 1, 3),
                (4, 1, 2, 3),
                (4, 0, 3, 3)):
            p = likelihoods[obs_to_idx[obs]]
            cumulative += p
            print(list(obs), p, cumulative, sep='\t')
        print()
        
        for variable in range(4):
            print('marginal distribution of variable %d:' % (variable + 1))
            distn = np.zeros(11)
            for obs in observations:
                p = likelihoods[obs_to_idx[tuple(obs)]]
                distn[obs[variable]] += p
            for i, m in enumerate(distn):
                print(i, m, sep='\t')
            mu = np.dot(distn, range(11))
            var = np.dot(distn, np.square(np.arange(11) - mu))
            std = np.sqrt(var)
            print('expected value:', mu)
            print('standard deviation:', std)
            print()
        print()


main()

output

  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
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
t: 10
-----

custom state probability requests:
[0, 0, 0, 10]	4.27269523375e-10
[6, 0, 0, 4]	0.00012185117912
[4, 0, 0, 6]	0.00876714664858

cumulative probabilities for constrained states:
[4, 0, 0, 6]	0.00876714664858	0.00876714664858
[4, 1, 0, 5]	0.00507950512166	0.0138466517702
[4, 0, 1, 5]	0.00247480222325	0.0163214539935
[4, 2, 0, 4]	0.00105475700432	0.0173762109978
[4, 1, 1, 4]	0.000188532804189	0.017564743802
[4, 0, 2, 4]	1.61715201033e-05	0.0175809153221
[4, 3, 0, 3]	5.01761227682e-06	0.0175859329344
[4, 2, 1, 3]	9.49455969192e-07	0.0175868823904
[4, 1, 2, 3]	8.93942169483e-08	0.0175869717846
[4, 0, 3, 3]	4.78150175753e-09	0.0175869765661

marginal distribution of variable 1:
0	0.649178233221
1	0.202352300698
2	0.0860071883895
3	0.0422680984811
4	0.0175869871342
5	0.00248435593013
6	0.000122636295128
7	1.99731186284e-07
8	1.19666830113e-10
9	3.0974428966e-14
10	2.92083903354e-18
expected value: 0.584677917955
standard deviation: 0.965852507559

marginal distribution of variable 2:
0	0.0210780436272
1	0.0458541453758
2	0.0913725137876
3	0.180224101344
4	0.288913531682
5	0.218097460737
6	0.153884858183
7	0.000574519641956
8	8.2508580153e-07
9	5.35302945887e-10
10	1.31687183721e-13
expected value: 3.94275029949
standard deviation: 1.45986075353

marginal distribution of variable 3:
0	0.621235915219
1	0.302342519445
2	0.0677440081271
3	0.00809823670744
4	0.00055876850802
5	2.05392116841e-05
6	1.27797453197e-08
7	1.20134503446e-12
8	3.94575617073e-17
9	5.11352182616e-22
10	2.20225778358e-27
expected value: 0.464463092599
standard deviation: 0.663272911013

marginal distribution of variable 4:
0	2.46865577089e-13
1	9.87776900735e-10
2	1.48275455143e-06
3	0.000990160757235
4	0.248815357757
5	0.495163529153
6	0.251180060624
7	0.00381961017802
8	2.96511191408e-05
9	1.46241484253e-07
10	4.27269523375e-10
expected value: 5.00810868995
standard deviation: 0.720729715281


t: 25
-----

custom state probability requests:
[0, 0, 0, 10]	7.85122364783e-09
[6, 0, 0, 4]	9.14250863738e-06
[4, 0, 0, 6]	0.00179032451313

cumulative probabilities for constrained states:
[4, 0, 0, 6]	0.00179032451313	0.00179032451313
[4, 1, 0, 5]	0.00103077852383	0.00282110303696
[4, 0, 1, 5]	0.000515185874101	0.00333628891107
[4, 2, 0, 4]	0.000230538717691	0.00356682762876
[4, 1, 1, 4]	5.66335223211e-05	0.00362346115108
[4, 0, 2, 4]	9.9984796116e-06	0.00363345963069
[4, 3, 0, 3]	2.06607964865e-06	0.00363552571034
[4, 2, 1, 3]	5.00586098348e-07	0.00363602629644
[4, 1, 2, 3]	8.53407967649e-08	0.00363611163723
[4, 0, 3, 3]	1.01334153622e-08	0.00363612177065

marginal distribution of variable 1:
0	0.732889748057
1	0.194209035319
2	0.0538331566312
3	0.0150871905072
4	0.00363613071316
5	0.000335483180635
6	9.24040123483e-06
7	1.51819829037e-08
8	9.28236580433e-12
9	2.45911305803e-15
10	2.3766165397e-19
expected value: 0.363414407614
standard deviation: 0.692932114012

marginal distribution of variable 2:
0	0.00976929611703
1	0.0276566005167
2	0.0747251236581
3	0.172732018089
4	0.292031097511
5	0.266305829669
6	0.155233763199
7	0.00154047475282
8	5.78679331099e-06
9	9.68809710466e-09
10	6.08870429492e-12
expected value: 4.13718872455
standard deviation: 1.32635888637

marginal distribution of variable 3:
0	0.628539203966
1	0.27267090024
2	0.0810796274669
3	0.0155777403443
4	0.00198978043154
5	0.000142406331116
6	3.41098120378e-07
7	1.21482679469e-10
8	1.4905396372e-14
9	7.11945270569e-19
10	1.1157389254e-23
expected value: 0.490236577028
standard deviation: 0.729565694382

marginal distribution of variable 4:
0	9.59283784417e-12
1	1.53515720275e-08
2	9.22086142183e-06
3	0.0024667990516
4	0.249533626069
5	0.491104858103
6	0.250220921167
7	0.00656052330496
8	0.000102880938165
9	1.14729208994e-06
10	7.85122364783e-09
expected value: 5.00916029081
standard deviation: 0.732671577687