How can I optimize coefficients of an arbitrary model? This might be terribly easy but I'm probably lacking the keywords to search for. Assume the following (dummy) data:
condition   session assumed_score  final_score
        1         1           110           95
        1         2           90            95
        1         3           80            95
        1         4           120           95
        1         5           130           95
        2         1           80            77
        2         2           90            77
        2         3           80            77
        2         4           50            77
        2         5           90            77
        3         1           60            35
        3         2           40            35
        3         3           30            35
        3         4           30            35
        3         5           40            35

In each testing condition, there are people taking part in 5 (or more, or less) sessions. They perform a task in each session. After all sessions, they give a final score (final_score) to their tasks. This is the ground truth. But we also (based on some measure) have assumed scores for each of the individual test sessions (assumed_score), based on observation of the test participants.
Now, I would like to be able—in the future—to predict the final score from the assumed scores. So, one "dumb" model would be to just average them:
$$avg_i = \frac{1}{N} \sum assumed_{ni}$$
where $N$ is the number of sessions and $i$ is the condition.
So now I have the averages, and I can build a simple linear model:
$$final = \beta_0 + \beta_1 avg$$
I could do this in R like so, and get a model that should fit the assumed scores to the final scores with least squares:
d = ddply(d, .(condition), transform, avg = mean(assumed_score))
lm(final_score ~ avg, data = d)

This works, because there are no coefficients in the average. 
But what if I wanted a more complex model? Say I'd like to give more weight on the sessions that have a lower score. Then my model would be:
$$final = \beta_0 + \left( \frac{1}{N} \sum_{n=1}^{N} assumed_n^{\beta} \right)^{\frac{1}{\beta}} $$
So here, I'd need to find an optimal parameter $\beta$ to fit my final scores against the assumed scores. But I have no idea how to do that with R.. I could try a fixed $\beta$, and then do:
custom <- function(d, beta = 0.1) {
  N = length(d)
  sum = 0
  for (a in d) {
    sum = sum + a^beta
  }
  avg = sum / N
  avg = avg^(1/beta)
}
d = ddply(d, .(condition), transform, custom = custom(assumed_score))
lm(final_score ~ custom, data = d) # => nice, this has lower R^2!

But that's hardly a good option to find a good value of $\beta$ without trying it manually.
How should I proceed?
 A: One problem is you haven't specified the error structure.
Your model isn't actually $final = \beta_0 + \left( \frac{1}{N} \sum_{n=1}^{N} assumed_n^{\beta} \right)^{\frac{1}{\beta}}$. The observed $y$-values won't equal the right hand side -- if that were the case you wouldn't need statistics. 
You might mean something like:
$final = \beta_0 + \left( \frac{1}{N} \sum_{n=1}^{N} assumed_n^{\beta} \right)^{\frac{1}{\beta}}+\epsilon$
where $\text{Var}(\epsilon)=\sigma^2I$.
In which case you probably want nonlinear least squares.
You would probably want to put your data in this (wide) form:
dw
  condition assumed_s1 assumed_s2 assumed_s3 assumed_s4 assumed_s5 final
1         1        110        120         90         90         30    95
2         2         90        130         80         60         30    77
3         3         80         80         50         40         40    35

since then
daffit <- nls(final ~ b0 +
   ((assumed_s1^b+assumed_s2^b+assumed_s3^b+assumed_s4^b+assumed_s5^b)/5)^(1/b),
   data=dw,start=list(b0=0,b=1))

> summary(daffit)

Formula: final ~ b0 + ((assumed_s1^b + assumed_s2^b + assumed_s3^b 
                          + assumed_s4^b + assumed_s5^b)/5)^(1/b)

Parameters:
   Estimate Std. Error t value Pr(>|t|)
b0  -17.232     31.944  -0.539    0.685
b     3.855     10.441   0.369    0.775

Residual standard error: 18.7 on 1 degrees of freedom

Number of iterations to convergence: 5 
Achieved convergence tolerance: 1.586e-06

If you know roughly what $\beta_0$ or $\beta$ should be (perhaps from some simpler or approximate model fit, for example, or from experience of typical values), use them as start values.

It's still possible to do it with a different number of sessions per condition ($N_i$ rather than $N$). 
With your particular function, you could still do it with nls since, for example, (if you keep $\beta>0\ $**) $0^\beta=0$, you can just add a column for $N_i$, pad out the number of terms to the maximum (with 0), then put $(\sum\{...\}/N_i)^{(1/\beta)}$ (where "$...$" has a sum of terms like now). 
** That can be done fairly simply via a reparameterization.
Some more complicated functions might require using more general optimization functions like optim, but you can get pretty clever with nls (cleverer than I have shown here).
Of course, none of this is really specific to R - many other stats packages have general optimization and/or flexible nonlinear least squares.
