1.4.1 Number of violations, \(f_k(SR_i, UR_i)\)

An important variable that shows up in the equations above is: \(f_k(SR_i, UR_i)\). This refers to the number of violations of constraint \(k\) incurred by the UR-SR pair, <\(UR_i, SR_i\)>. For example, the UR-SR pair \(<\)/ur-2/, [sr2-2]\(>\) picks up two violations of Con2, so \(f_{Con2}(\text{sr2-2}, \text{ur-2}) = 2\) (yellow box).

1.4.2 Penalty

Next, we move on to the penalty. Penalty has also been called “harmony” in other descriptions of the MaxEnt model. To calculate the penalty for the UR-SR pair <\(UR_i, SR_i\)>, we take the weighted sum of its constraint violations: \[\text{penalty for <}UR_i, SR_i\text{> } = \text{ } \Sigma_{k=1}^m w_k f_k(SR_i, UR_i)\] (Recall that \(w_k\) is the weight of constraint \(k\), and \(m\) is the number of constraints.)

Let’s apply this equation to calculate the penalty for the UR-SR pair \(<\)/ur-2/, [sr2-2]\(>\). This pair picks up one violation of Con1 and two violations of Con2. The weighted sum (green box) is thus:

\[\begin{align} \text{penalty for </ur-2/, [sr2-2]>} &= \Sigma_{k=1}^2 w_k f_k(\text{sr2-2, ur-2}) \\ &= weight_{Con1} \times 1 + weight_{Con2} \times 2 \\ &= 1.6 \times 1 + 0.8 \times 2 \\ &= 3.2 \end{align}\]

The penalties for other UR-SR pairs are similarly calculated. For example, the penalty for the UR-SR pair \(<\)/ur-2/, [sr2-1]\(>\) is:

\[\begin{align} \text{penalty for </ur-2/, [sr2-1]>} &= \Sigma_{k=1}^2 w_k f_k(\text{sr2-1, ur-2}) \\ &= weight_{Con1} \times 0 + weight_{Con2} \times 1 \\ &= 1.6 \times 0 + 0.8 \times 1 \\ &= 0.8 \end{align}\]

Notice that \(m=2\) because there are two constraints.

1.4.3 MaxEnt score

The next column gives the MaxEnt score for each UR-SR pair. \[\text{MaxEnt score} = exp(-\text{penalty})\] The Reader may have also encountered this term being called the “potential” in other descriptions of MaxEnt. Look at the MaxEnt score column of Figure 2 to see this equation applied.

1.4.4 Normalization, \(Z\)

We then need to normalize the MaxEnt score in order to get probabilities. The next column provides the \(Z\) value for each UR. \(Z\) is short for Zustandssumme, “sum over states”, and has also been known as the normalization constant, or the canonical partition function. \(Z\) is the sum of the MaxEnt scores of all UR-SR pairs arising from a particular UR.

Here’s a simplified way to write the equation for the normalization constant of \(UR_i\): \[Z(UR_i) = \Sigma_{j=1}^l m_j\] where \(l\) is the number of UR-SR pairs whose UR is \(UR_i\), and \(m_j\) is the MaxEnt score of the \(j^{th}\) UR-SR pair whose UR is \(UR_i\).

Let’s apply this equation to calculate the normalization constant for /ur-2/ (blue box):

\[\begin{align} Z(\text{ur-2}) &= \Sigma_{i=1}^2m_i \\ &= 0.499 + 0.041 \\ &= 0.490 \end{align}\]

Here, \(l = 2\) because there are only two UR-SR pairs whose UR is /ur-2/. We basically add the MaxEnt scores of these two pairs to get \(Z(\text{ur-2})\).

1.4.5 Predicted probability

The predicted probability is shown in the final column. The probability of a particular SR for a given UR, \(UR_i\), is this pair’s MaxEnt score divided by the normalization constant associated with \(UR_i\):

\[P(SR_i|UR_i) = \frac{\text{MaxEnt score of <}UR_i, SR_i\text{>}}{Z(UR_i)}\]

Let’s apply this equation to calculate the probability of [sr2-2] given /ur-2/ (magenta box):

\[\begin{align} P(\text{sr2-2}|\text{ur-2}) &= \frac{\text{MaxEnt score of <ur-2, sr2-2>}}{Z(\text{ur-2})} \\ &= \frac{0.449}{0.490} \\ &= 0.083 \end{align}\]

Congratulations! You now know how to calculate the probabilities of an SR given its UR when provided with constraint weights. Next, we’ll discuss some characteristics of the probability distribution and of the penalty score.

2.3.1 Defining log likelihood

When no priors are used, maximizing the objective function \(J(w)\) is equivalent to maximizing the log likelihood of the weights \(LL(w|D)\) alone. Since we’ve not used any priors so far, in all the preceding cases, the sole aim of the learner was to seek weights that maximized the log likelihood alone.

The log likelihood of the weights is:

\[LL(w|D) = \Sigma_{i=1}^n ln P(y_i|x_i; w)\]

Let’s unpack the terms of this equation. On the left-hand side, we have \(LL(w|D)\), where

• \(LL(w|D)\) is the log likelihood of the weights \(w\) given the observed data \(D\)
• \(w\) is the vector of constraint weights
• \(D\) is the training data

On the right-hand side:

• \(P(y_i|x_i; w)\) is the probability of realizing underlying \(x_i\) as surface \(y_i\) under the weights \(w\)
• \(n\) is the number of observed data points
  • i.e. the sum of the frequency column in the input
• \(y_i\) is the observed surface form of the \(i^{th}\) data point, and \(x_i\) is the input form that \(y_i\) arose from

Notice how the weights of each individual constraint (which are found in \(w\)) are the only parameters available for the model to adjust.

2.3.2 What guides the weight changes?: Maximizing log likelihood

Recall that the learner takes in a distribution over training examples, \(D\). For the data set we’ve been using, \(D\) is: \(<\)Input1, Output1-1\(>\): 1, \(<\)Input1, Output1-2\(>\): 1, \(<\)Input2, Output2-1\(>\): 1. \(<\)Input2, Output2-2\(>\) is not part of \(D\) because it isn’t observed.

In the absence of priors, learner’s sole goal is to maximize the log likelihood alone. In other words, the learner seeks the weights that produce a distribution over these three training examples that matches the observed frequencies as closely as possible. The learner adjusts the weights with this goal in mind.

Table 1: Log likelihood of M1’s & M2’s weights.

Let’s imagine we have two models, M1 & M2, whose weights (whatever they may be) produce the predicted distributions P(output|input) shown in the final two columns of Table 1. For ease of comparison, the observed relative frequencies are shown in the third-to-last column. Right away, we can see that M1 matches the observed frequencies much better than M2 does. In other words, the observed data are more probable under the weights of M1 than they are under the weights of M2. Though we haven’t yet applied the equation above to calculate the log likelihood of these two model’s weights, we can expect that M1’s weights will produce a higher log likelihood than M2’s.

Now we’ll apply the log likelihood equation to calculate the log likelihood of the weights of M1 & M2. Recall that log likelihoods may range from negative infinity to 0, with higher log likelihoods indicating that a model’s weights produce a better fit to the observed data. Here’s the log likelihood equation again:

\[LL(w|D) = \Sigma_{i=1}^n ln P(y_i|x_i; w)\]

Applying the equation to M1 (\(w_{M1}\) refers to M1’s weights): \[\begin{align} LL(w_{M1}|D) &= \Sigma_{i=1}^3 ln P(y_i|x_i; w) \\ &= ln(.5) + ln(.5) + ln(1) \\ &= -1.386 \end{align}\]

Applying the equation to M2 (\(w_{M2}\) refers to M2’s weights): \[\begin{align} LL(w_{M2}|D) &= \Sigma_{i=1}^3 ln P(y_i|x_i; w) \\ &= ln(.25) + ln(.75) + ln(1) \\ &= -2.367 \end{align}\]

M1’s weights produce a higher log likelihood than M2’s, which indicates that M1 is a better model than M2. This is because M1’s weights produce a predicted distribution over observed UR-SR pairs that is closer to the observed frequencies than M2’s weights do. (Notice that only observed data are included in the calculations above. Despite Output2-2 getting 50% of the predicted probability in M2, it does not contribute to the log likelihood.)

Hence, aiming to maximize the probability of the training data results in finding weights that produce a predicted distribution that is as close to the observed frequencies as possible. For example, we can imagine the learner starting with the weights that produce a distribution that is very distant from the frequencies of the observed training data. With its goal of finding weights that maximize the log likelihood, the learner moves to the weights of M2, and eventually arrives at the weights of M1, where training is halted and M1’s weights are the weights of the final trained model.

3.3.1 Ceiling log likelihood

In §2.2 we noted that the log likelihood attained by simple_model was in fact at ceiling.

# Get log likelihood of simple model (no biases)
simple_model$loglik
#> [1] -1.386296

Each data set has a theoretical maximum log likelihood, which can be calculated. Recall that the data set had 3 data points (i.e. observed outputs). (Frequency shows up in the third column.)

# View the data that was used to train `simple_model`
read.table(data_file, header=FALSE, sep='\t')
#>                                               V1
#> 1 Input,Output,Frequency,Constraint1,Constraint2
#> 2                         Input1,Output1-1,1,1,0
#> 3                         Input1,Output1-2,1,0,1
#> 4                         Input2,Output2-1,1,0,0
#> 5                         Input2,Output2-2,0,0,1

• For Input1, there were 2 outputs, each of which surfaced once.
  • Observed \(P(\text{Output1-1}|\text{Input1}) = 0.5\)
  • Observed \(P(\text{Output1-2}|\text{Input1}) = 0.5\)
• For Input2, 1 output surfaced once.
  • Observed \(P(\text{Output2-1}|\text{Input2}) = 1\)

Let’s calculate the theoretical ceiling log likelihood. Here’s the equation for the log likelihood, repeated for easy reference:

\[LL(w|D) = \Sigma_{i=1}^n ln P(y_i|x_i; w)\]

A trained model that matches the observed frequencies perfectly would have attained the ceiling log likelihood. We’ll imagine just such a perfect trained model that matches the observed frequencies perfectly, and calculate its log likelihood: \[\begin{align} LL(w|D) &= \Sigma_{i=1}^3 ln P(y_i|x_i; w)\\ &= ln(.5) + ln(.5) + ln(1) \\ &= -1.386 \end{align}\]

Since the log likelihood attained by our perfect trained model matches that of simple_model, we conclude that simple_model also attained the ceiling log likelihood. In other words, the distribution produced by simple_model’s weights matches the observed frequencies of the training data perfectly.

3.3.2 Regularization as a strategy to combat overfitting

Trained models that fit the training data too well often cannot generalize to unseen data. Such models are said to overfit the training data. One strategy to reduce overfitting is to include a regularization term. The intended effect of regularization is to keep the weights as low as possible.

To put this into practice, we’ll make use of the optional arguments mu & sigma of the function optimize_weights(). To keep weights low, we’ll set mu = 0 & sigma = 1. This means that we believe that every constraint should:

• Have an ideal weight = 0
• Have standard deviation = 1

Since every constraint has the same ideal weight & the same standard deviation (tolerance for weight discrepancy), we can make use of the optional arguments mu & sigma to input \(\mu\) and \(\sigma\) just once for all constraints.

# Train regularized model with mu=0 & sigma=1
regularized_model <- optimize_weights(
  data_frame, mu=0, sigma=1
)

Inspecting the trained regularized model, we now see a value for the attribute bias_params:

# Take a took at the regularized model trained with scalar biases
regularized_model
#> $name
#> [1] "data_frame"
#> 
#> $weights
#> Constraint1 Constraint2 
#>   0.1051951   0.3163665 
#> 
#> $loglik
#> [1] -1.944845
#> 
#> $k
#> [1] 2
#> 
#> $n
#> [1] 3
#> 
#> $bias_params
#>     Constraint    Mu Sigma
#>         <char> <num> <num>
#> 1: Constraint1     0     1
#> 2: Constraint2     0     1

We can retrieve this attribute alone:

# View the bias parameters we used during model training
regularized_model$bias_params
#>     Constraint    Mu Sigma
#>         <char> <num> <num>
#> 1: Constraint1     0     1
#> 2: Constraint2     0     1

Constraints are ordered in rows. The first column gives the \(\mu\) and the second column the \(\sigma\) for the constraint in question.

Let’s take a look at the trained weights for the regularized model:

# Get learned weights for the regularized model trained with scalar biases
regularized_model$weights
#> Constraint1 Constraint2 
#>   0.1051951   0.3163665

The trained regularized model has learned that Constraint 1 has weight = 0.105, and Constraint 2 has weight = 0.316. Notice that the weights for both constraints are now closer to 0 (cf. the unbiased model where trained weights for both constraints were 14.2). We attained this by including a prior term that tagged on a price to pay for weights whose value did not exactly match the ideal weight. In particular the price to be paid increased the further a constraint’s weight was from the ideal weight \(0\).

Notice also that the log likelihood of the regularized model’s weights is below the ceiling of \(-1.386\) (\(-1.945 < -1.386\)):

# Get value of objective function for trained regularized model
regularized_model$loglik
#> [1] -1.944845

The regularized model provides a worse fit to the observed data than the unregularized model does; but the regularized model is probably able to better generalize to unseen data.

Retrieving the weight of 1 constraint

Let’s imagine that we’re only interested in a tracking the weight of a subset of the constraints. For example, we’re interested only in Constraint 1, and want to track the weight learned for this constraint in the unbiased and the simple bias (regularized) model. We use index 1 to pick out the first constraint from the named list your_model_name$weights. Similarly, if we were interested in the second constraint, we’d use index 2, and so on.

# Get weight of Constraint 1 (simple model)
cons1_noBias <- simple_model$weights[1]

# Get weight of Constraint 1 (regularized model)
cons1_simpleBias <- regularized_model$weights[1]

# Compare learned weights of Constraint 1 in the unbiased model...
# ...and in the regularized one
print(paste("In the unbiased model, Constraint 1's  weight is:",
            sprintf(cons1_noBias, fmt = '%#.3f')))
#> [1] "In the unbiased model, Constraint 1's  weight is: 13.521"
print(paste("In the regularized model, Constraint 1's  weight is:",
            sprintf(cons1_simpleBias, fmt = '%#.3f')))
#> [1] "In the regularized model, Constraint 1's  weight is: 0.105"

3.4.1 The phenomenon: Productivity of saltation (White 2017)

Consider the following artificial grammar learning experiment (White 2017). Participants were trained on either a p→v alternation (experimental condition) or a b→v alternation (control condition). During the response phase, participants were tested on two alternations. In the first, they were tested on producing the alternation they were trained on (“test”). In the second, they were tested on extending the alternation to the other unseen pair (“extension”). That is, for those trained on the p→v alternation, the b→v alternation was considered the “untrained alternation/extension”, and vice versa.

Experimental result to be modeled:

Extension to untrained sounds occurred more frequently in the experimental condition (trained on p→v; extended to b→v) at 73% than in the control condition (trained on b→v; extended to p→v) at 20%. The percentage alternation by participants in the response phase is shown below.

#>           Expt cdn Control cdn
#> trained         95          90
#> untrained       73          20

Given that the perceptual distance between [b] & [v] is smaller than that between [p] & [v], substantive effects arising from perceptual distance appear to play an important role.

3.4.2 Modeling the extension of p→v (trained) to b→v (untrained)

So how can we introduce these substantive factors into the model? Intuitively, we want to make it more difficult to map [p] to [v] than it is to map [b] to [v]. We will use *Map(x,y) constraints (Zuraw 2007, 2013), which prohibit the mapping between sound x and sound y. We give *Map(p,v) a higher preferred weight than *Map(b,v), which means that the mapping from [p] to [v] is more strongly prohibited than the mapping from [b] to [v]. Following White, we’ll set \(\mu\)’s of these constraints to the following preferred weights:

• \(\mu\) of *Map(b,v) = 1.3
• \(\mu\) of *Map(p,v) = 3.65

Markedness constraints *V[-cont]V and *V[-voi]V drive the alternation from stops to [v] intervocally. These have a preferred weight of 0 to reflect that intervocalic stops and voiceless sounds were allowed in the participants’ native languages.

The training data for the experimental condition includes only the p→v alternation, reflecting the training phase encountered by the subjects in the experimental condition:

# Get path to data file (experimental condition)
white_salt_data_file <- system.file(
  "extdata", "sample_whiteSalt_data_file.txt", package = "maxent.ot"
)

# File is in OTSoft format, so first convert it to data frame
white_salt_data_frame <- otsoft_tableaux_to_df(white_salt_data_file)

# View training data (experimental condition)
white_salt_data_frame
#>   Input Output Frequency X.V..cont.V X.V..voi.V X.Map.b.v. X.Map.p.v.
#> 1   VpV    VpV         0           1          1          0          0
#> 2   VpV    VbV         0           1          0          0          0
#> 3   VpV    VvV        18           0          0          0          1

We now train the bias model for the experimental condition, making use of the argument mu. We want the value of \(\mu\) to be 0, 0, 1.3 and 3.65 for the four ordered constraints [*V[-cont]V, *V[-voi]V, *Map(b,v), *Map(p,v)] respectively. We have to combine these four values into a vector using the c() function (i.e. c(0, 0, 1.3, 3.65)). Here, we choose to have the same tolerance for deviation from the preferred weight for all constraints, so we’ll use the argument sigma.

# Fit model with preferred weights for experimental condition
white_salt_bias_model = optimize_weights(
  white_salt_data_frame, mu=c(0, 0, 1.3, 3.65), sigma=sqrt(0.6)
)

# View trained weights (expt cdn with bias)
white_salt_bias_model$weights
#> X.V..cont.V  X.V..voi.V  X.Map.b.v.  X.Map.p.v. 
#>   2.5877273   0.8014478   1.3000000   1.0622343

Now let’s train an unbiased model for the same experimental condition. We’ll average out the preferred weights for the two *Map constraints, which results in each of these constraint having preferred weight = 2.27.

# Fit unbiased model for experimental condition
white_salt_noBias_model = optimize_weights(
  white_salt_data_frame, 
  mu=c(0, 0, 2.27, 2.27), sigma=sqrt(0.6)
)

# View trained weights (expt cdn, no bias)
white_salt_noBias_model$weights
#> X.V..cont.V  X.V..voi.V  X.Map.b.v.  X.Map.p.v. 
#>   2.0607951   0.6888778   2.2700000   0.2091823

At this point, we’ve essentially trained the model on only the p→v alternation. Next, we’ll see how well our trained models perform on both the test task (p→v) as well as on the extension task (b→v). Notice how both these phases (training and response) reflect the experimental setup with human subjects.

We’ll need a new test file that includes the observed frequencies for both the test (p→v) and extension (b→v) alternations. We’ll pass both the models trained above as well as the new test file to the function predict_probabilities(). We’ll take a closer look at both the function predict_probabilities(), and the structure of the test file when we introduce predict_probabilities() in §4 below. For now, we’ll just focus our attention on the predicted probability calculated by this function.

# Get path to test file (experimental condition)
white_salt_test_file <- system.file(
  "extdata", "sample_whiteSalt_test_file.txt", package = "maxent.ot"
)

# File is in OTSoft format, so convert first
white_salt_test_df <- otsoft_tableaux_to_df(white_salt_data_file)

# View test data (experimental condition)
white_salt_test_df
#>   Input Output Frequency X.V..cont.V X.V..voi.V X.Map.b.v. X.Map.p.v.
#> 1   VpV    VpV         0           1          1          0          0
#> 2   VpV    VbV         0           1          0          0          0
#> 3   VpV    VvV        18           0          0          0          1

# Predict probabilities with weights trained with bias (expt cdn)
predict_probabilities(
  white_salt_test_df, white_salt_bias_model$weights
)
#> $loglik
#> [1] -4.930533
#> 
#> $predictions
#>   Input Output Freq X.V..cont.V X.V..voi.V X.Map.b.v. X.Map.p.v.  Predicted
#> 1   VpV    VpV    0           1          1          0          0 0.07420978
#> 2   VpV    VbV    0           1          0          0          0 0.16539619
#> 3   VpV    VvV   18           0          0          0          1 0.76039403
#>   Observed       Error
#> 1        0  0.07420978
#> 2        0  0.16539619
#> 3        1 -0.23960597

# Predict probabilities with weights trained without bias (expt cdn)
predict_probabilities(
  white_salt_test_df, white_salt_noBias_model$weights
)
#> $loglik
#> [1] -3.811101
#> 
#> $predictions
#>   Input Output Freq X.V..cont.V X.V..voi.V X.Map.b.v. X.Map.p.v. Predicted
#> 1   VpV    VpV    0           1          1          0          0 0.0637862
#> 2   VpV    VbV    0           1          0          0          0 0.1270289
#> 3   VpV    VvV   18           0          0          0          1 0.8091849
#>   Observed      Error
#> 1        0  0.0637862
#> 2        0  0.1270289
#> 3        1 -0.1908151

For ease of reference, the predicted percentage of alternation for the trained biased and unbiased models are summarized in the table below:

#> [1] "Predicted % of alternation during response phase (experimental condition)"
#>                 Biased model Unbiased model
#> Trained (95%)             91             93
#> Untrained (73%)           78             45

Recall that during the response phase, participants performed the following two alternations at the following rates:

• [p] → [v] (trained “test”): 95%
• [b] → [v] (untrained “extension”): 73%

Overall, the biased model provided a better fit to the observed data. Predicted probability differed by no more than 5.5%, with the greatest discrepancy showing up for the untrained [b] → [v] alternation.

The unbiased model, in contrast, performed poorly. In particular, the untrained [b] → [v] alternation appeared to be close to chance (45%).

3.4.3 Modeling the extension of b→v (trained) to p→v (untrained)

Let’s now train a biased and an unbiased model for the control condition. The training data includes only the b→v alternation, reflecting the training phase of the control condition:

# Get path to data file (control condition)
white_control_data_file <- system.file(
  "extdata", "sample_whiteCtrl_data_file.txt", package = "maxent.ot"
)

white_control_df <- otsoft_tableaux_to_df(white_control_data_file)

# View training data (control condition)
white_control_df
#>   Input Output Frequency X.V..cont.V X.V..voi.V X.Map.b.v. X.Map.p.v.
#> 1   VbV    VpV         0           1          1          0          0
#> 2   VbV    VbV         0           1          0          0          0
#> 3   VbV    VvV        18           0          0          1          0

# Fit model with preferred weights for control condition
white_control_bias_model = optimize_weights(
  white_control_df,
  mu=c(0, 0, 1.3, 3.65), sigma=sqrt(0.6)
)

# View trained weights (control cdn with bias)
white_control_bias_model$weights
#> X.V..cont.V  X.V..voi.V  X.Map.b.v.  X.Map.p.v. 
#>   1.9372214   0.6600614   0.0000000   3.6500000

# Fit unbiased model for control condition
white_control_noBias_model = optimize_weights(
  white_control_df,
  mu=c(0, 0, 2.27, 2.27), sigma=sqrt(0.6)
)

# View trained weights (control cdn, no bias)
white_control_noBias_model$weights
#> X.V..cont.V  X.V..voi.V  X.Map.b.v.  X.Map.p.v. 
#>   2.0607951   0.6888778   0.2091823   2.2700000

As with the experimental condition above, we’ll now subject the two trained models to the response phase, which includes the test (b→v) and extension (p→v) tasks.

# Get path to test file (control condition)
white_control_test_file <- system.file(
  "extdata", "sample_whiteCtrl_test_file.txt", package = "maxent.ot"
)

# File is in OTSoft format, so convert first
white_control_test_df <- otsoft_tableaux_to_df(white_control_test_file)

# View test data (experimental condition)
white_control_test_df
#>   Input Output Frequency X.V..cont.V X.V..voi.V X.Map.b.v. X.Map.p.v.
#> 1   VpV    VpV       0.8           1          1          0          0
#> 2   VpV    VvV       0.2           0          0          0          1
#> 3   VbV    VbV       0.1           1          0          0          0
#> 4   VbV    VvV       0.9           0          0          1          0

# View test data (control condition)
read.table(white_control_test_file, header=FALSE, sep='\t')
#>    V1  V2  V3         V4        V5        V6        V7
#> 1          NA *V[-cont]V *V[-voi]V *Map(b-v) *Map(p-v)
#> 2          NA        Mk1       Mk2       Mp1       Mp2
#> 3 VpV VpV 0.8          1         1         0         0
#> 4     VvV 0.2          0         0         0         1
#> 5 VbV VbV 0.1          1         0         0         0
#> 6     VvV 0.9          0         0         1         0

# Predict probabilities with weights trained with bias (control cdn)
predict_probabilities(
  white_control_test_df, white_control_bias_model$weights
)
#> $loglik
#> [1] -0.838242
#> 
#> $predictions
#>   Input Output Freq X.V..cont.V X.V..voi.V X.Map.b.v. X.Map.p.v. Predicted
#> 1   VpV    VpV  0.8           1          1          0          0 0.7412963
#> 2   VpV    VvV  0.2           0          0          0          1 0.2587037
#> 3   VbV    VbV  0.1           1          0          0          0 0.1259534
#> 4   VbV    VvV  0.9           0          0          1          0 0.8740466
#>   Observed       Error
#> 1      0.8 -0.05870366
#> 2      0.2  0.05870366
#> 3      0.1  0.02595343
#> 4      0.9 -0.02595343

# Predict probabilities with weights trained with bias (control cdn)
predict_probabilities(
  white_control_test_df, white_control_noBias_model$weights
)
#> $loglik
#> [1] -1.196516
#> 
#> $predictions
#>   Input Output Freq X.V..cont.V X.V..voi.V X.Map.b.v. X.Map.p.v. Predicted
#> 1   VpV    VpV  0.8           1          1          0          0 0.3823294
#> 2   VpV    VvV  0.2           0          0          0          1 0.6176706
#> 3   VbV    VbV  0.1           1          0          0          0 0.1356836
#> 4   VbV    VvV  0.9           0          0          1          0 0.8643164
#>   Observed       Error
#> 1      0.8 -0.41767061
#> 2      0.2  0.41767061
#> 3      0.1  0.03568365
#> 4      0.9 -0.03568365

For ease of reference, the predicted percentage of alternation for the trained biased and unbiased models are presented in the table below:

#> [1] "Predicted % of alternation during response phase (control condition)"
#>                 Biased model Unbiased model
#> Trained (90%)             87             87
#> Untrained (20%)           26             61

During the testing phase, participants performed the following two alternations at the following rates:

• [b] → [v] (trained “test”): 90%
• [p] → [v] (untrained “extension”): 20%

The biased model matches the observed data relatively well. The unbiased model, in contrast, greatly overgeneralizes to the untrained p→v alternation.

3.4.4 Conclusion

Overall, the biased model, which incorporated substantive biases, performed much better than the unbiased one (no substantive influences) at modeling human subject responses during the response phase. The unbiased model both greatly underpredicted extension to the b→v alternation in the experimental condition and greatly overpredicted the extension to the p→v alternation in the control condition. In contrast, the biased model neither vastly under- nor overpredicted on any of the response tasks.

For the biased model, we introduced different preferred weights for separate constraints by utilizing the mu argument of the optimize_weights() function. These preferred weights were chosen to represent the substantive effects that arose from perceptual distances. The resulting trained biased model was able to match the observed frequencies of alternations produced by human subjects during the response phase much more closely than its unbiased counterpart in all experimental conditions, thus demonstrating the utility of incorporating substantive biases for linguistic modeling.

5.3.1 Hungarian backness harmony

In Hungarian, the dative suffix has two allomorphs: [-nɔk] and [-nɛk]. A front rounded vowel anywhere in the stem triggers the front allomorph. When the closest vowel is a back vowel, the back allomorph is triggered. However, if a back vowel is separated from the suffix by one or more neutral vowels, variation is observed. Likewise, if a stem consists of only neutral vowels, variation is observed.

In addition, several unnatural data patterns were observed. For example, stems with a final bilabial stop favored the front allomorph. Such patterns were deemed unnatural because (non-dorsal) consonants are generally thought not to have any effect on vowel harmony.

Constraints associated with vowel-triggered harmony were deemed “natural” while those associated with consonant-triggered ones were deemed “unnatural”. In addition, monosyllabic Ci:C stems take a back allomorph (historically, this was a back vowel), so the constraint USE_BACK/Ci:C__ was also considered unnatural.

HZSL were interested in assessing whether the natural and unnatural constraints governing the zones of variation were under- or over-learned. They designed a wug test where Hungarian speakers were presented with a wug stem as well as the allomorphs [-nɔk] and [-nɛk]. The speakers had to pick one of the two suffixes to go with the stem, so this was a 2AFC task.

5.3.2 Challenges with assessing under-/over-learning

When one wishes to assess whether humans under- or over-learn certain constraints, it is natural to want to compare the lexicon-trained weights (Grammar L) directly to the wug-test-trained weights (Grammar W). However, such direct comparisons can be problematic for several reasons:

1. The lexicon and the wug test are based on different types of words

• Wug tests are often carefully designed to capture a only a particular part of the grammar that the linguist is interested in
• In other words, entire types of words may be missing from the wug test.
• For example, HZSL found that in the lexicon, all stems whose final nucleus contained a back vowel (B-type stems) always picked a back suffix.
  • Lexicon model (Grammar L): Agree(back, local) has very high weight.
• However, such stems were very uninteresting to HZSL, who were interested in obtaining wug results for stems in the zones of variation.
• HZSL excluded such B-type stems from their wug test:
  • Wug model (Grammar W): Agree(back, local) has extremely low weight, since the learner never encountered any candidate that violated this constraint.

2. Even if all word types are included in the wug test, the proportions between the different types must be identical to that of the lexicon in order for the wug & lexicon grammars to be comparable.
3. Human subject responses are depolarized if experimental task is 2AFC

5.3.3 A rigorous approach to assessing under-/over-learning

Instead of directly comparing a constraint’s weight in Grammar L to its weight in Grammar W, here’s what HZSL did:

1. HSZL simulated a group of idealized Hungarian speakers who had learned the weights from the lexicon perfectly (i.e. they interact with the world using Grammar L).
2. This group of idealized Hungarian speakers then performed the very same wug test that the group of human subjects did, randomly picking the back or front suffix for each stem according to their grammar (i.e. the weights in Grammar L).
3. The simulated wug test result was then fit to produce a grammar, Grammar S.

• Notice how this mirrored the human wug test results that were fit to produce the grammar, Grammar W.
• In other words, both Grammar S and Grammar W were trained on data sets that had the exact same types and proportion of stems.

Comparing Grammar S to Grammar W:

• We know that Grammar S arose from the wug test results of a group of idealized speakers using Grammar L, because HZSL made these simulated speakers perform the wug test using Grammar L.
• We also have Grammar W, whose weights we know arose from the wug test results of human Hungarian speakers, but we don’t know what grammar these speakers used to perform the wug test.
• If indeed the human Hungarian speakers had learned their lexicon perfectly, then they should have used Grammar L to perform the wug test.
• Consequently, Grammar S and Grammar W should not differ.
• However, if a particular constraint was under-learned by the human speakers, then its weight should be lower in Grammar W than it is in Grammar S.
• The same applied for over-learned constraints, which would have a higher weight in Grammar W than in Grammar S.

A huge collection of Grammar S:

• Since the stem types in the wug test were those in the zones of variation, that lone wug test performed by one group of idealized speakers might produce very different suffix choices than the one performed by the human speakers due to chance alone.
• To mitigate this, HZSL ran 10,000 groups of simulated idealized speakers.
• For each constraint C, they now had a distribution of its expected idealized weight (from 10,000 Grammar S’s)
• This approach, estimating the value of a variable via random sampling, is called the Monte Carlo method.
• If constraint C’s weight from Grammar W fell in the left tail of the idealized C distribution, they interpreted this to mean that C was under-learned.
• If constraint C’s weight from Grammar W fell in the right tail of the idealized C distribution, then C was over-learned.

Let’s take stock of where we are. In §5.3.2, we saw that there were three challenges to directly comparing the lexicon and wug grammars. Simulating idealized speakers performing the same wug task as human speakers solved the first two issues. The third issue, simulating depolarized 2AFC wug responses, will be taken up in the next section (§5.3.4), where the hyperparameter temperature will be used to effect depolarization.

5.3.4 Generating depolarized predictions for a 2AFC task

How would a simulated speaker go about picking a suffix for a wug stem? For each stem, the simulated speaker’s grammar (i.e. Grammar L) generates a probability distribution over the two suffix options. This probability distribution guides the simulated speaker’s suffix choice. For example, if the probability distribution is:

• back suffix: 90%
• front suffix: 10%

then 9 out of 10 times, the simulated speaker will pick the back suffix.

We will utilize the predict_probabilities() function to generate the probability distribution over the two suffix options of a given stem:

• test_file: Hungarian wug test data
  • Since we want to get the predicted probabilities over the suffix options of the various wug stems.
• constraint_weights: Grammar L
  • Since these idealized speakers have learned the lexicon perfectly, we’ll use the weights that were learned by training on the lexicon.
• temparature: T=1.5
  • To capture the depolarization effect for 2AFC tasks, we need \(T>1\)
  • Following HZSL, we’ll specifically use \(T=1.5\). See HZSL (2009) for details on why they chose this value.

# Get paths to sample Hungarian wug data
hu_file <- system.file(
  "extdata", "sample_hu_wug.txt", package = "maxent.ot"
)

hu_data <- otsoft_tableaux_to_df(hu_file)

# We'll use the weights learned by training on the lexicon as reported by HZSL on (p.848).
lex_wts_ls <- list(5.39, 1.69, 1.46, 3.00, 4.04, 2.45, 1.08, .91, 1.73, 2.42)

# Make a named list (optional)
names(lex_wts_ls) <- list(
  'AGREE(back,nonlocal)',   'AGREE(front,local)',
  'AGREE(nonhigh_front,local)', 'AGREE(low_front,local)',
  'AGREE(double_front,local)',  'USE_FRONT/bilabial__',
  'USE_FRONT/[+cor,+son]__', 'USE_FRONT/sibilant__',
  'USE_FRONT/CC__', 'USE_BACK/[C0i:C0]__'
)

# Convert to double object (required type as argument of predict_probabilities())
lex_wts <- as.double(lex_wts_ls)

# Predict probabilities with temperature set to 1.5
# Notice also that "UTF-8" is an option for encoding, which comes in useful here
hu_pred <- predict_probabilities(
  hu_data, lex_wts, temperature = 1.5, encoding = "UTF-8"
)

Let’s take a look at some of the predicted wug probabilities produced by Grammar L:

# Let's view some of the predicted probabilities at T=1.5
# We've dropped the constraint violations for easier viewing
# Both predicted and observed probability are conditioned over URs
# e.g. Observed probability in rows 3 & 4 don't sum to 1
# because the stem [étt] appears twice in this data set.
head(hu_pred$predictions[, -4:-(ncol(hu_pred$predictions)-3)])
#>      Input      Output Freq  Predicted Observed       Error
#> 1 Szólingy Szólingynak    1 0.78807041      1.0 -0.21192959
#> 2 Szólingy Szólingynek    0 0.21192959      0.0  0.21192959
#> 3      étt      éttnak    0 0.01860365      0.0  0.01860365
#> 4      étt      éttnek    1 0.48139635      0.5 -0.01860365
#> 5 Csazenip Csazenipnak    0 0.13470305      0.0  0.13470305
#> 6 Csazenip Csazenipnek    1 0.86529695      1.0 -0.13470305

# For curiosity's sake, let's compare:
# With default temperature=1...
# ... predicted probability is more polarized
hu_pred_t1 <- predict_probabilities(
  hu_data, lex_wts, encoding = "UTF-8"
)
head(hu_pred_t1$predictions[, -4:-(ncol(hu_pred_t1$predictions)-3)])
#>      Input      Output Freq   Predicted Observed        Error
#> 1 Szólingy Szólingynak    1 0.877611113      1.0 -0.122388887
#> 2 Szólingy Szólingynek    0 0.122388887      0.0  0.122388887
#> 3      étt      éttnak    0 0.003769867      0.0  0.003769867
#> 4      étt      éttnek    1 0.496230133      0.5 -0.003769867
#> 5 Csazenip Csazenipnak    0 0.057866955      0.0  0.057866955
#> 6 Csazenip Csazenipnek    1 0.942133045      1.0 -0.057866955

5.3.5 Monte-carlo simulations based on the predicted probabilities above

We’ve used Grammar L to generate predicted probabilities for wug suffixes at \(T=1.5\), and stored the probability distributions in the object hu_pred. Next, we use these probability distributions to guide our simulated subjects as they pick suffixes for wug stems.

We use the function monte_carlo_weights() to simulate 500 groups of idealized Hungarian speakers picking suffixes for the wug stems. This function only requires two arguments: the prediction data frame (we use hu_pred), and the number of simulations.

This function produces a complete wug test response by randomly choosing one suffix option for each trial in hu_pred based on the the predicted probabilities in hu_pred. (This function automatically converts predicted P(SR|UR) in hu_pred to predicted P(SR|trial).) It then calls optimize_weights() to fit this simulated wug response (i.e. to produce Grammar S), and records the weights learned. It repeats this process for the specified number of simulations. It returns a data frame with \(m\) rows and \(n\) columns, where \(m\) is the number of simulations and \(n\) is the number of constraints. That is, each row represents an instance of Grammar S.

The following code was used to produce the 500 simulated complete wug test responses and train their corresponding grammars (Grammar S’s):

# Learn weights for 500 simulated wug tests
# monte_carlo_weights(hu_pred, 500)

Each complete wug test consisted of 1703 trials, so the process took a very long time to run. To save time, we saved the 500 sets of constraint weights (500 Grammar S’s) and present them below:

# Get path to the saved constraint weights trained on the 500 simulated wug tests
hu_500simul_wts_path <- system.file(
  "extdata", "hu_500simuls_wts.txt", package = "maxent.ot"
)

# Store constraint weights in object `hu_500simul_wts`
hu_500simul_wts <- read.table(hu_500simul_wts_path, header=TRUE, sep='\t')
hu_500simul_wts <- data.frame(hu_500simul_wts)

# Rename columns
names(hu_500simul_wts) <- c(
  'AGREE(back,nonlocal)',   'AGREE(front,local)',
  'AGREE(nonhigh_front,local)', 'AGREE(low_front,local)',
  'AGREE(double_front,local)',  'USE_FRONT/bilabial__',
  'USE_FRONT/[+cor,+son]__', 'USE_FRONT/sibilant__',
  'USE_FRONT/CC__', 'USE_BACK/[C0i:C0]__'
)

# View the first 6 rows of hu_500simul_wts
# i.e. the first 6 Grammar S's
head(hu_500simul_wts)
#>   AGREE(back,nonlocal) AGREE(front,local) AGREE(nonhigh_front,local)
#> 1             3.383119          1.0876934                  0.9360368
#> 2             3.444765          0.9890333                  0.8731836
#> 3             3.588559          1.1637950                  0.8753390
#> 4             3.323187          1.1766615                  0.8087822
#> 5             3.509050          0.9975359                  0.9905781
#> 6             3.394041          0.5826657                  1.2360608
#>   AGREE(low_front,local) AGREE(double_front,local) USE_FRONT/bilabial__
#> 1               1.838278                  2.624190             1.593205
#> 2               2.107615                  2.506750             1.703648
#> 3               2.180713                  3.195825             1.874271
#> 4               2.085677                  2.435280             1.824118
#> 5               1.919777                  2.912883             1.605194
#> 6               1.818039                  2.976685             1.635315
#>   USE_FRONT/[+cor,+son]__ USE_FRONT/sibilant__ USE_FRONT/CC__
#> 1               0.6393167            0.7599629       0.966830
#> 2               0.6472897            0.8210698       1.216772
#> 3               0.2137989            0.4346787       1.244254
#> 4               0.3060865            0.6178177       1.099641
#> 5               0.7514726            0.7526335       1.074058
#> 6               1.1377657            0.8029462       1.302554
#>   USE_BACK/[C0i:C0]__
#> 1           1.4764315
#> 2           2.1874098
#> 3           0.8430519
#> 4           2.2069378
#> 5           1.8524475
#> 6           1.6132052

Some housekeeping

Notice that duplicate UR entries are present in the prediction data frame hu_pred. This is because the set-up of the monte_carlo_weights() function crucially depends upon predicted P(SR|trial) in order to randomly give one SR in each trial a frequency of 1. Consecutive rows are considered to belong to the same trial if they have the same UR, but not a repeated SR. When creating your prediction data frame, it is important to ensure that trials are not collapsed into a single entry. The caution against duplicate UR entries highlighted in §4.1 doesn’t apply here because the Error column of the hu_pred file isn’t used by monte_carlo_weights().

The process took a long time to run so we performed only 500 (rather than 10,000) simulations. The Reader is invited to uncomment the code below if they would like to try running the monte_carlo_weights() function.

# Learn weights for 5 simulated wug tests
#hu_simul_wts <- monte_carlo_weights(hu_pred, 5)
#hu_simul_wts

5.3.6 Are unnatural constraints under-learned?

Now that we have a distribution of idealized constraint weights for each of the 10 constraints (from 500 Grammar S’s), let’s find out whether the weights trained on human wug responses (Grammar W) were under-, over-, or perfectly learned.

First, we’ll have to fit a model to the human wug responses (i.e. produce Grammar W):

# Fit model for human wug response
human_wug_model <- optimize_weights(hu_data)

Let’s view the constraint weights of Grammar W in relation to the average idealized weights from the 500 Grammar S’s:

# Weights trained on human wug response
human_wt <- human_wug_model$weights

# Average idealized wug weights
ideal_wt <- colMeans(hu_500simul_wts)

# Create data frame
wt_df <- data.frame(human_wt, ideal_wt)

# View weights trained on human wug responses & average idealized wug weights
print(wt_df, digits = 3)
#>                            human_wt ideal_wt
#> AGREE.back.nonlocal.          3.581    3.635
#> AGREE.front.local.            2.223    1.155
#> AGREE.nonhigh_front.local.    0.950    0.985
#> AGREE.low_front.local.        1.052    2.007
#> AGREE.double_front.local.     1.939    2.724
#> USE_FRONT.bilabial__          1.037    1.631
#> USE_FRONT...cor..son.__       0.433    0.725
#> USE_FRONT.sibilant__          0.368    0.605
#> USE_FRONT.CC__                0.688    1.163
#> USE_BACK..C0i.C0.__           0.805    1.616

All unnatural constraints (those beginning with USE…) appear to be under-learned by humans. For example, USE_FRONT/bilabial__ has a lower weight when trained on human wug responses (1.04) than the average of the weights obtained by training on simulated idealized wug responses (1.63).

The plots below visualize where the weights trained on human wug tests (blue vertical dashed line) fall in relation to the distribution of the weights trained on the simulated idealized wug tests (histogram).

# Plot weights for the five unnatural constraints by looping through them
# For i-th constraint

for (i in 6:10) {
  hist(hu_500simul_wts[[i]], breaks=15, xlab="weight", main = colnames(hu_500simul_wts[i]))
  abline(v=human_wug_model$weights[i], col="red", lwd=3, lty=2)
}

6.1.1 Likelihood ratio test: Background

Two requirements must be met in order to perform the likelihood ratio test. First, exactly two models must be compared. Second, these models must be nested.

Here, we have two data sets small_data & large_data (thus meeting the first criterion). large_data has the very same constraints (C1 & C3) that small_data has, plus an additional constraint, C2. Thus, the nested criterion is met.

# Get paths to toy data files
# This file has two constraints
small_file <- system.file(
  "extdata", "sample_data_file_small.txt", package = "maxent.ot"
)
small_df <- otsoft_tableaux_to_df(small_file)

# This file has three constraints
large_file <- system.file(
  "extdata", "sample_data_file_large_otsoft.txt", package = "maxent.ot"
)
large_df <- otsoft_tableaux_to_df(large_file)

# View small_data
small_df
#>    Input    Output Frequency Constraint1 Constraint3
#> 1 Input1 Output1-1         1           1           1
#> 2 Input1 Output1-2         1           0           0
#> 3 Input2 Output2-1         1           0           1
#> 4 Input2 Output2-2         0           0           0

# View large_data
large_df
#>    Input    Output Frequency Constraint1 Constraint2 Constraint3
#> 1 Input1 Output1-1         1           1           0           1
#> 2 Input1 Output1-2         1           0           1           0
#> 3 Input2 Output2-1         1           0           0           1
#> 4 Input2 Output2-2         0           0           1           0

We’ll train both models and store the trained models in the objects small_model & large_model:

# Fit weights to both data sets
small_model <- optimize_weights(small_df)
large_model <- optimize_weights(large_df)

Now we’ll perform the likelihood ratio test on the two trained models. To use the lrt method of compare_models(), we’ll need three arguments:

1. object storing results of first model
2. object storing results of second model
3. method = 'lrt'

The models can appear in either order, but they must both appear before the method argument.

# Compare models using the likelihood ratio test
compare_models(small_model, large_model, method='lrt')
#>         description   chi_sq k_delta   p_value
#> 1 large_df~small_df 1.386292       1 0.2390322

This method returns a data frame with a single row and the following columns:

• description: The names of the two models being compared, with the larger model appearing first.
• chi_sq: chi-squared value.
• k_delta: degrees of freedom (difference in the number of constraints between the two model).
• p-value: p-value.

The p-value was \(0.239\). This indicates that the difference in model fit to data was not significant at the \(.05\) significance level. Hence, we are not justified in including the additional constraint, C2. The model with fewer constraints is the better one.

6.1.2 Likelihood ratio test: Explained

The likelihood ratio, \(LR\), is calculated as follows:

\[LR = 2(LL_2 - LL_1)\]

where \(LL_2\) is the log likelihood of the larger model (the one with more constraints) and \(LL_1\) is the log likelihood of the smaller model.

The likelihood ratio \(LR\) is set as the \(X^2\)-value, which is used to perform a \(X^2\)-test. The \(X^2\)-value should range between 0 & positive infinity, with a larger value indicating a greater improvement in model fit. The difference in the number of constraints is set as the degrees of freedom. The p-value tells us whether the difference in likelihood between the two models is significant. If it turns out to be significant, then we are justified in including the additional constraints in order to achieve the better model fit.

This test runs on the assumption that the larger model indeed produces a better fit to the data. If this criterion isn’t met, there is no way in which the larger model is better than the smaller one since

1. the larger model does worse on model fit, and
2. even with more constraints, the larger model couldn’t fit the data better than the smaller model

Mathematically, this translates to having a negative \(X^2\)-value, which is uninterpretable. An error message will be produced and no \(X^2\)-test will be performed if such a scenario is encountered.

6.1.3 Likelihood ratio won’t work for non-nested models

We’ll attempt to perform a likelihood ratio test on two non-nested models. We re-use small_data, which has constraints C1 & C3. We introduce a new data set, large_data_b, which has constraints C1, C2 & C4.

# Get paths to toy data files
# This file has three constraints
large_file_b <- system.file(
  "extdata", "sample_data_file_large_b.txt", package = "maxent.ot"
)
large_df_b <- otsoft_tableaux_to_df(large_file_b)

# View small_data
small_df
#>    Input    Output Frequency Constraint1 Constraint3
#> 1 Input1 Output1-1         1           1           1
#> 2 Input1 Output1-2         1           0           0
#> 3 Input2 Output2-1         1           0           1
#> 4 Input2 Output2-2         0           0           0

# View large_data_b
large_df_b
#>    Input    Output Frequency Constraint1 Constraint2 Constraint4
#> 1 Input1 Output1-1         1           1           0           0
#> 2 Input1 Output1-2         1           0           1           0
#> 3 Input2 Output2-1         1           0           0           1
#> 4 Input2 Output2-2         2           0           1           0

# Fit weights to data set large_data_b 
large_model_b <- optimize_weights(large_df_b)

The program is unable to perform the likelihood ratio test, and will produce an error message.

# Compare models using the likelihood ratio test
compare_models(small_model, large_model_b, method='lrt')
#> Error in compare_models(small_model, large_model_b, method = "lrt"): Constraint names in the smaller model are not a subset of those in the larger model. Models must be nested to apply the likelihood ratio test.

6.2.1 AIC: Background

The AIC value for a particular model is calculated as follows:

\[AIC = 2k - 2LL\] where \(k\) is the number of constraints, and \(LL\) is the log likelihood of the current model.

A smaller AIC score indicates a better model. Both the number of constraints and the log likelihood contribute to the AIC score. All else being equal, having a smaller number of constraints or having a higher log likelihood is indicative of a better model. Recall that log likelihood (\(-\infty < LL \leq 0\)) quantifies model fit.

6.2.2 Compare models using the AIC method

Two or more models may be passed to compare_models() when using the aic method. To use the aic method of compare_models(), we’ll need \(m+1\) arguments, where \(m\) is the number of models that we want compared.

• First argument: object storing results of first model
• Second argument: object storing results of second model
• …
• \(m^{th}\) argument: object storing results of \(m^{th}\) model
• \((m+1)^{th}\) argument: method = 'aic'

The models can appear in any order, but they must all appear before the method argument.

The returned data frame contains the same number of rows as models passed in. The models are sorted in increasing order of AIC (i.e. best first). The data frame has the following columns:

• model: The name of the model.
• k: The number of constraints.
• aic: The model’s AIC score.
• aic.delta: The difference between the current model’s AIC score and the AIC score of the “best” model.
• aic.wt: The model’s AIC weight: This reflects the relative likelihood (or conditional probability) that this model is the “best” model in the set.
• cum.wt: The cumulative sum of AIC weights up to and including this model.
• ll: The log likelihood of this model.

# Compare models using AIC
compare_models(small_model, large_model, method='aic')
#>      model k      aic aic.delta    aic.wt    cum.wt        ll
#> 1 small_df 2 8.158883 0.0000000 0.5761171 0.5761171 -2.079442
#> 2 large_df 3 8.772591 0.6137077 0.4238829 1.0000000 -1.386295

In the data frame above, the small model has an AIC score of 8.16 while the large model has an AIC score of 8.77. This indicates that the small model is better.

How much better is the small model compared to the large model? The column aic.wt gives the probability that the current model is the best one of the entire set of models passed in. So the small model has a 57.6% chance of being the best model, and the large model has a 43.7% chance of being the best model.

6.2.3 AIC weight (aic.wt) explained

The AIC weight of a model is calculated by applying the softmax function to the negative halved-distance between current model’s AIC score and the AIC score of the best model:

\[W_i = \frac{exp(-0.5\delta_i)}{\Sigma_{j=1}^m exp(-0.5\delta_j)}\] where \(\delta_i\) is the difference in AIC score between model \(i\) and the best model, and \(m\) is the number of models being compared.

We’ll break this equation down further, and walk through its application step-by-step. The values of the intermediate terms we’ll use (e.g. \(\delta\); the raw AIC weight, \(r\)) are shown in Figure 5.

Fig 5: Calculating AIC weights: Intermediate steps

6.3.1 BIC: Background

The BIC is an alternative measure of assessing the trade-off between the number of constraints and model fit. The BIC replaces the AIC’s first term (\(2k\)) with \(k \ln(n)\).

\[BIC = k \ln(n) - 2 LL\]

where \(k\) is the number of constraints, \(n\) is the sample size, and \(LL\) is the log likelihood of the current model.

As with the AIC, a smaller BIC score indicates a better model.

For both the AIC and the BIC, the first term penalizes models that have more constraints, \(k\). However, the number of constraints, \(k\), is scaled differently in these two measures. For the BIC, \(k\) is scaled by \(\ln(n)\), where \(n\) is the number of samples. As the number of observed tokens increases, the penalty accrued by models that have more constraints increases (albeit at a decreasing rate since \(k\) is scaled by \(\ln(n)\) rather than \(n\)). In contrast, the AIC scales \(k\) with the constant \(2\), so sample size has no effect on the penalty accrued by having more constraints, \(k\). At \(7\) observed tokens and below, a model’s BIC score is smaller than its AIC score. At \(8\) observed tokens and above, a model’s BIC score is larger than its AIC score (since \(\ln(7) < 2 < \ln(8)\)).

A caveat: Though unlikely to occur, it is advisable to check that the sample size is greater than \(1\). At a sample size of \(1\), the effect of the number of constraints, \(k\), on the BIC score is canceled out. The BIC score becomes directly proportional to the log likelihood alone (\(BIC = k \ln(1) -2LL = -2LL\)).

6.3.2 Compare models using the BIC method

As with the AIC, two or more models may be passed to compare_models() when using the bic method. The arguments passed in to compare_models for the bic method are identical to that for the aic method, except that the method argument should now be: method='bic'.

Hence, to use the bic method of compare_models(), we’ll need \(m+1\) arguments, where \(m\) is the number of models that we want compared.

• First argument: object storing results of first model
• Second argument: object storing results of second model
• …
• \(m^{th}\) argument: object storing results of \(m^{th}\) model
• \((m+1)^{th}\) argument: method = 'bic'

The models can appear in any order, but they must all appear before the method argument.

# Compare models using BIC
compare_models(small_model, large_model, method='bic')
#>      model k n      bic bic.delta    bic.wt    cum.wt        ll
#> 1 large_df 3 3 6.068428   0.00000 0.5358981 0.5358981 -1.386295
#> 2 small_df 2 3 6.356108   0.28768 0.4641019 1.0000000 -2.079442

The returned data frame from the BIC method has the same structure as that of the AIC, with one additional column:

• n: The number of samples in the data that the model is fit to.

The bic and bic.wt columns are the BIC counterparts to the aic and aic.wt columns in the AIC data frame.

Looking at the data frame above, we see that the small model (smaller BIC score) again performs better than the large model (large BIC score). The conditional probability that the small model is the best of all models passed to compare_models() is 53.6%.

6.4.1 AIC-C: Background

The Akaike Information Criterion for small sample sizes (AIC-C) yet another alternative method of comparing the trade-off between the number of constraints and the model fit. Like the BIC, AIC-C also takes the sample size into account.

The AIC-C modifies the AIC by introducing a third term.

\[AIC_c = 2k - 2 LL + \frac{2k^2 + 2k}{n-k-1}\]

As before, \(k\) is the number of constraints, \(n\) is the sample size, and \(LL\) is the log likelihood of the current model. This newly introduced third term further increases the cost of having additional constraints. As \(n\) approaches infinity, the third term’s denominator approaches 0, so the AIC-C becomes equivalent to the AIC for large sample sizes.

6.4.2 Compare models using the AIC-C method

As with the AIC and the BIC, two or more models may be passed to compare_models() when using the aic_c method. The arguments passed in to compare_models for the aic_c method are identical to that for the aic and bic methods, except that the method argument should now be: method='aic_c'.

Hence, to use the bic method of compare_models(), we’ll need \(m+1\) arguments, where \(m\) is the number of models that we want compared.

• First argument: object storing results of first model
• Second argument: object storing results of second model
• …
• \(m^{th}\) argument: object storing results of \(m^{th}\) model
• \((m+1)^{th}\) argument: method = 'aic_c'

The models can appear in any order, but they must all appear before the method argument.

# Get paths to toy data files
# This file has 2 constraints
small_file_c <- system.file(
  "extdata", "sample_data_file_small_c.txt", package = "maxent.ot"
)
small_df_c <- otsoft_tableaux_to_df(small_file_c)

# This file has 3 constraints
large_file_c <- system.file(
  "extdata", "sample_data_file_large_c.txt", package = "maxent.ot"
)
large_df_c <- otsoft_tableaux_to_df(large_file_c)

# Fit weights to both data sets 
small_model_c <- optimize_weights(small_df_c) 
large_model_c <- optimize_weights(large_df_c)

# Compare models using AIC-C
compare_models(small_model_c, large_model_c, method='aic_c')
#>        model k n     aicc aicc.delta    aicc.wt    cum.wt        ll
#> 1 small_df_c 2 6 16.31777   0.000000 0.97375561 0.9737556 -4.158883
#> 2 large_df_c 3 6 23.54518   7.227416 0.02624439 1.0000000 -2.772591

The returned data frame from the AIC-C method has the same structure as that of the BIC. The aicc and aicc.wt columns are the AIC-C counterparts to the aic and aic.wt columns in the AIC data frame respectively.

In the data frame above, the small model (smaller AIC-C score) performs better than the large model (larger AIC-C score). The conditional probability that the small model is the best of all the models that were passed to compare_models() is \(97.4%\).

6.4.3 A potential pitfall

Here’s the formula for AIC-C again:

\[AIC_c = 2k - 2 LL + \frac{2k^2 + 2k}{n-k-1}\]

When \(n-k \leq 1\) for a particular model, the third term’s denominator is 0. This results in an infinitely large AIC-C value. In the example below, this issue crops up for both the large and the small models.

# Problematic data set: all models have n-k <= 1

# Compare models using AIC-C
compare_models(small_model, large_model, method='aic_c')
#>      model k n aicc aicc.delta aicc.wt cum.wt        ll
#> 1 small_df 2 3  Inf        NaN     NaN    NaN -2.079442
#> 2 large_df 3 3  Inf        NaN     NaN    NaN -1.386295

In the data frame above, none of the models produce a valid AIC-C score. Hence the probability that one particular model is the best cannot be calculated (aicc.wt is NaN “Not a Number”).

One way to side-step this issue is to scale the sample size up. For example, the problematic data sets had their sample sizes increased by a factor of 2 to create the non-problematic data sets:

# View problematic data set
# sample size: n = 3
# maximum number of constraints for valid aic_c score: 3-2 = 1
# actual maximum number of constraints used: 3
large_df
#>    Input    Output Frequency Constraint1 Constraint2 Constraint3
#> 1 Input1 Output1-1         1           1           0           1
#> 2 Input1 Output1-2         1           0           1           0
#> 3 Input2 Output2-1         1           0           0           1
#> 4 Input2 Output2-2         0           0           1           0

# View non-problematic data set
# sample size: n = 6
# maximum number of constraints for valid aic_c score: 6-2 = 4
# actual maximum number of constraints used: 3
large_df_c
#>    Input    Output Frequency Constraint1 Constraint2 Constraint3
#> 1 Input1 Output1-1         2           1           0           1
#> 2 Input1 Output1-2         2           0           1           0
#> 3 Input2 Output2-1         2           0           0           1
#> 4 Input2 Output2-2         0           0           1           0

In the case where only one (i.e. large_model_d) of multiple models produces an invalid AIC-C score, the model with the invalid AIC-C score will have an aicc.wt of 0. The aicc.wt for the other models with valid AIC-C scores will be calculated as per normal (see data frame below).

# 1 Problematic data set: large_data_d has n-k <= 1
# 2 Non-problematic data sets

# sample size: n = 5
# maximum number of constraints for valid aic_c score: 5-2 = 3

# Get paths to toy data files
# This file has 2 constraints: Good
small_file_d <- system.file(
  "extdata", "sample_data_file_small_d.txt", package = "maxent.ot"
)
small_df_d <- otsoft_tableaux_to_df(small_file_d)

# This file has 3 constraints: Good
med_file_d <- system.file(
  "extdata", "sample_data_file_med_d.txt", package = "maxent.ot"
)
med_df_d <- otsoft_tableaux_to_df(med_file_d)

# This file has 4 constraints: Problematic
large_file_d <- system.file(
  "extdata", "sample_data_file_large_d.txt", package = "maxent.ot"
)
large_df_d <- otsoft_tableaux_to_df(large_file_d)

# Fit weights for all data sets 
small_model_d <- optimize_weights(small_df_d)
med_model_d <- optimize_weights(med_df_d)
large_model_d <- optimize_weights(large_df_d)

# Compare models using AIC-C
compare_models(small_model_d, med_model_d, large_model_d, method='aic_c')
#>        model k n     aicc aicc.delta      aicc.wt    cum.wt        ll
#> 1 small_df_d 2 5 16.73012    0.00000 9.999513e-01 0.9999513 -3.365058
#> 2   med_df_d 3 5 36.59167   19.86156 4.865154e-05 1.0000000 -3.295837
#> 3 large_df_d 4 5      Inf        Inf 0.000000e+00 1.0000000 -3.295837