diff --git a/cmd/lb/main.go b/cmd/lb/main.go
index 7d24271..24a9c29 100644
--- a/cmd/lb/main.go
+++ b/cmd/lb/main.go
@@ -216,7 +216,7 @@ func main() {
 			os.Exit(1)
 		}
 		var weights, reps []float64
-		var date []time.Time
+		var dates []time.Time
 		isbw := review.IsBodyweight(ps[0].Exercise, ps[0].Type)
 
 		for _, p := range (ps) {
@@ -227,16 +227,14 @@ func main() {
 				if p.RIR[i] == 0 {
 					weights = append(weights, w)
 					reps = append(reps, float64(p.Reps[i]))
-					date = append(date, p.Times[i].UTC())
+					dates = append(dates, p.Times[i].UTC())
 					fmt.Printf("Set: weight %0.0f, reps %0.0f\n", w, float64(p.Reps[i]))
 				}
 			}
 		}
 		//a, b := review.FitPowerLaw(weights, reps, date, 14.0)
-		est, err := review.Fit(weights, reps, date,
-				review.WithAllModels(),
-				review.WithHalfLife(14.0),
-			)
+		est := review.NewEstimator(review.WithHalfLife(14.0))
+		err = est.Fit(weights, reps, dates)
 		if err != nil {
 			fmt.Printf("Error estimating performance: %v\n", err)
 			os.Exit(1)
@@ -246,7 +244,7 @@ func main() {
 			if isbw {
 				adj = ps[len(ps)-1].Bodyweight
 			}
-			fmt.Printf("%d: %0.0f\n", i, est.EstimateMaxWeight(float64(i)) - adj)
+			fmt.Printf("%d: %0.0f\n", i, est.EstimateMaxWeight(time.Now(), float64(i)) - adj)
 		}
 	case "predict":
 		if flag.NArg() != 3 {
diff --git a/review/estimate.go b/review/estimate.go
index f2d028c..a255c08 100644
--- a/review/estimate.go
+++ b/review/estimate.go
@@ -2,199 +2,170 @@ package review
 
 import (
 	"errors"
-	"fmt"
 	"math"
+	"sort"
 	"time"
 
 	"gonum.org/v1/gonum/optimize"
 )
 
-// Estimator encapsulates a fitted model and exposes estimation methods.
+// Estimator is the main interface for incremental, time-aware strength estimation.
 type Estimator interface {
-	Estimate1RM() float64
-	EstimateReps(targetWeight float64) float64
-	EstimateMaxWeight(nReps float64) float64
-	ModelType() string
-	Params() []float64
+	Fit(weight, reps []float64, dates []time.Time) error
+	Estimate1RM(t time.Time) float64
+	EstimateReps(t time.Time, targetWeight float64) float64
+	EstimateMaxWeight(t time.Time, nReps float64) float64
+	Params(t time.Time) []float64
 }
 
-// Supported model types
-const (
-	ModelPowerLaw   = "powerlaw"
-	ModelLinear     = "linear"
-	ModelExponential = "exponential"
-)
+// --- Functional Options ---
 
-// FitOption is a functional option for configuring the Fit process.
-type FitOption func(*fitConfig)
-
-// fitConfig holds configuration for fitting.
-type fitConfig struct {
-	modelTypes   []string
-	halfLifeDays float64
+type estimatorConfig struct {
+	modelType   string
+	halfLife    float64
+	smoothAlpha float64 // Exponential smoothing factor (0 < alpha <= 1)
 }
 
-// WithModel specifies which model(s) to fit. If multiple, Fit selects the best.
-func WithModel(models ...string) FitOption {
-	return func(cfg *fitConfig) {
-		cfg.modelTypes = models
+type EstimatorOption func(*estimatorConfig)
+
+// WithModel sets the model type (currently only "powerlaw" is implemented).
+func WithModel(model string) EstimatorOption {
+	return func(cfg *estimatorConfig) {
+		cfg.modelType = model
 	}
 }
 
-// WithAllModels configures Fit to try all built-in model types.
-func WithAllModels() FitOption {
-	return func(cfg *fitConfig) {
-		cfg.modelTypes = []string{ModelPowerLaw, ModelLinear, ModelExponential}
+// WithHalfLife sets the half-life for time weighting in curve fitting.
+func WithHalfLife(days float64) EstimatorOption {
+	return func(cfg *estimatorConfig) {
+		cfg.halfLife = days
 	}
 }
 
-// WithHalfLife sets the half-life (in days) for time weighting.
-func WithHalfLife(days float64) FitOption {
-	return func(cfg *fitConfig) {
-		cfg.halfLifeDays = days
+// WithSmoothingAlpha sets the exponential smoothing factor for parameter smoothing.
+func WithSmoothingAlpha(alpha float64) EstimatorOption {
+	return func(cfg *estimatorConfig) {
+		cfg.smoothAlpha = alpha
 	}
 }
 
-// WithModelSelection is an alias for WithModel, for clarity.
-func WithModelSelection(models []string) FitOption {
-	return WithModel(models...)
+// --- Estimator Implementation ---
+
+type timePoint struct {
+	date time.Time
+	a    float64
+	b    float64
 }
 
-// Default settings
-const defaultHalfLife = 30.0
-var defaultModelTypes = []string{ModelPowerLaw}
+type estimatorImpl struct {
+	cfg        estimatorConfig
+	data       []timePoint // sorted by date
+	smoothedA  []float64   // smoothed a for each timePoint
+	smoothedB  []float64   // smoothed b for each timePoint
+}
 
-// Fit fits the specified model(s) to the data and returns an Estimator.
-// If multiple models are specified, Fit selects the best based on residual sum of squares.
-func Fit(weight, reps []float64, dates []time.Time, opts ...FitOption) (Estimator, error) {
-	if len(weight) != len(reps) || len(weight) != len(dates) {
-		return nil, errors.New("weight, reps, and dates must have the same length")
-	}
-	if len(weight) < 2 {
-		return nil, errors.New("at least two data points are required")
-	}
-
-	// Apply options
-	cfg := &fitConfig{
-		modelTypes:   defaultModelTypes,
-		halfLifeDays: defaultHalfLife,
+// NewEstimator creates a new Estimator with the given options.
+func NewEstimator(opts ...EstimatorOption) Estimator {
+	cfg := estimatorConfig{
+		modelType:   "powerlaw",
+		halfLife:    30.0,
+		smoothAlpha: 0.3,
 	}
 	for _, opt := range opts {
-		opt(cfg)
+		opt(&cfg)
 	}
-	if len(cfg.modelTypes) == 0 {
-		cfg.modelTypes = defaultModelTypes
+	return &estimatorImpl{
+		cfg: cfg,
 	}
-
-	// Fit each model and select the best (lowest residual)
-	var best Estimator
-	var bestResidual float64 = math.Inf(1)
-	now := time.Now()
-
-	for _, model := range cfg.modelTypes {
-		var est Estimator
-		var residual float64
-		var err error
-
-		switch model {
-		case ModelPowerLaw:
-			est, residual, err = fitPowerLaw(weight, reps, dates, now, cfg.halfLifeDays)
-		case ModelLinear:
-			est, residual, err = fitLinear(weight, reps, dates, now, cfg.halfLifeDays)
-		case ModelExponential:
-			est, residual, err = fitExponential(weight, reps, dates, now, cfg.halfLifeDays)
-		default:
-			return nil, fmt.Errorf("unknown model type: %s", model)
-		}
-		if err != nil {
-			continue // Skip models that fail to fit
-		}
-		if residual < bestResidual {
-			best = est
-			bestResidual = residual
-		}
-	}
-	if best == nil {
-		return nil, errors.New("no model could be fitted to the data")
-	}
-	return best, nil
 }
 
-// --- Model Implementations ---
+// Fit adds new data and updates the parameter time series and smoothing.
+func (e *estimatorImpl) Fit(weight, reps []float64, dates []time.Time) error {
+	if len(weight) != len(reps) || len(weight) != len(dates) {
+		return errors.New("weight, reps, and dates must have the same length")
+	}
+	// Add new data points
+	for i := range weight {
+		e.data = append(e.data, timePoint{
+			date: dates[i],
+			a:    math.NaN(), // to be filled in
+			b:    math.NaN(),
+		})
+	}
+	// Sort all data points by date
+	sort.Slice(e.data, func(i, j int) bool {
+		return e.data[i].date.Before(e.data[j].date)
+	})
 
-// PowerLawEstimator: w = a * reps^b
-type PowerLawEstimator struct {
-	a, b        float64
-	halfLife    float64
-	modelType   string
-	residualSum float64
+	// For each time point, fit the model to all data up to that point
+	for i := range e.data {
+		var w, r []float64
+		var d []time.Time
+		for j := 0; j <= i; j++ {
+			w = append(w, weight[j])
+			r = append(r, reps[j])
+			d = append(d, dates[j])
+		}
+		a, b := fitPowerLaw(w, r, d, e.cfg.halfLife)
+		e.data[i].a = a
+		e.data[i].b = b
+	}
+
+	// Smooth the parameter time series
+	e.smoothedA = exponentialSmoothing(extractA(e.data), e.cfg.smoothAlpha)
+	e.smoothedB = exponentialSmoothing(extractB(e.data), e.cfg.smoothAlpha)
+	return nil
 }
 
-func (e *PowerLawEstimator) Estimate1RM() float64 {
-	return e.a * math.Pow(1, e.b)
+// Estimate1RM returns the smoothed 1RM estimate at time t.
+func (e *estimatorImpl) Estimate1RM(t time.Time) float64 {
+	a, b := e.smoothedParamsAt(t)
+	return a * math.Pow(1, b)
 }
-func (e *PowerLawEstimator) EstimateReps(targetWeight float64) float64 {
-	if e.a == 0 || e.b == 0 {
+
+// EstimateReps returns the predicted number of reps at a given weight and time.
+func (e *estimatorImpl) EstimateReps(t time.Time, targetWeight float64) float64 {
+	a, b := e.smoothedParamsAt(t)
+	if a == 0 || b == 0 {
 		return 0
 	}
-	return math.Pow(targetWeight/e.a, 1/e.b)
-}
-func (e *PowerLawEstimator) EstimateMaxWeight(nReps float64) float64 {
-	return e.a * math.Pow(nReps, e.b)
-}
-func (e *PowerLawEstimator) ModelType() string { return e.modelType }
-func (e *PowerLawEstimator) Params() []float64 { return []float64{e.a, e.b} }
-
-// LinearEstimator: w = a + b*reps
-type LinearEstimator struct {
-	a, b        float64
-	halfLife    float64
-	modelType   string
-	residualSum float64
+	return math.Pow(targetWeight/a, 1/b)
 }
 
-func (e *LinearEstimator) Estimate1RM() float64 {
-	return e.a + e.b*1
+// EstimateMaxWeight returns the predicted max weight for a given number of reps at time t.
+func (e *estimatorImpl) EstimateMaxWeight(t time.Time, nReps float64) float64 {
+	a, b := e.smoothedParamsAt(t)
+	return a * math.Pow(nReps, b)
 }
-func (e *LinearEstimator) EstimateReps(targetWeight float64) float64 {
-	if e.b == 0 {
-		return 0
+
+// Params returns the smoothed model parameters at time t.
+func (e *estimatorImpl) Params(t time.Time) []float64 {
+	a, b := e.smoothedParamsAt(t)
+	return []float64{a, b}
+}
+
+// --- Internal Helpers ---
+
+// smoothedParamsAt returns the smoothed parameters for the closest time point <= t.
+func (e *estimatorImpl) smoothedParamsAt(t time.Time) (float64, float64) {
+	if len(e.data) == 0 {
+		return 0, 0
 	}
-	return (targetWeight - e.a) / e.b
-}
-func (e *LinearEstimator) EstimateMaxWeight(nReps float64) float64 {
-	return e.a + e.b*nReps
-}
-func (e *LinearEstimator) ModelType() string { return e.modelType }
-func (e *LinearEstimator) Params() []float64 { return []float64{e.a, e.b} }
-
-// ExponentialEstimator: w = a * exp(b * reps)
-type ExponentialEstimator struct {
-	a, b        float64
-	halfLife    float64
-	modelType   string
-	residualSum float64
-}
-
-func (e *ExponentialEstimator) Estimate1RM() float64 {
-	return e.a * math.Exp(e.b*1)
-}
-func (e *ExponentialEstimator) EstimateReps(targetWeight float64) float64 {
-	if e.a == 0 || e.b == 0 {
-		return 0
+	idx := sort.Search(len(e.data), func(i int) bool {
+		return !e.data[i].date.Before(t)
+	})
+	if idx == 0 {
+		return e.smoothedA[0], e.smoothedB[0]
 	}
-	return math.Log(targetWeight/e.a) / e.b
+	if idx >= len(e.data) {
+		return e.smoothedA[len(e.data)-1], e.smoothedB[len(e.data)-1]
+	}
+	return e.smoothedA[idx-1], e.smoothedB[idx-1]
 }
-func (e *ExponentialEstimator) EstimateMaxWeight(nReps float64) float64 {
-	return e.a * math.Exp(e.b*nReps)
-}
-func (e *ExponentialEstimator) ModelType() string { return e.modelType }
-func (e *ExponentialEstimator) Params() []float64 { return []float64{e.a, e.b} }
 
-// --- Fitting Functions ---
-
-// fitPowerLaw fits w = a * reps^b
-func fitPowerLaw(weight, reps []float64, dates []time.Time, now time.Time, halfLifeDays float64) (Estimator, float64, error) {
+// fitPowerLaw fits a power law model to the data.
+func fitPowerLaw(weight, reps []float64, dates []time.Time, halfLifeDays float64) (a, b float64) {
+	now := dates[len(dates)-1]
 	params := []float64{max(weight), -0.1}
 	problem := optimize.Problem{
 		Func: func(x []float64) float64 {
@@ -203,16 +174,9 @@ func fitPowerLaw(weight, reps []float64, dates []time.Time, now time.Time, halfL
 	}
 	result, err := optimize.Minimize(problem, params, nil, nil)
 	if err != nil {
-		return nil, 0, err
+		return 0, 0
 	}
-	residual := weightedResidualsPowerLaw(result.X, weight, reps, dates, now, halfLifeDays)
-	return &PowerLawEstimator{
-		a:          result.X[0],
-		b:          result.X[1],
-		halfLife:   halfLifeDays,
-		modelType:  ModelPowerLaw,
-		residualSum: residual,
-	}, residual, nil
+	return result.X[0], result.X[1]
 }
 
 func weightedResidualsPowerLaw(params, weight, reps []float64, dates []time.Time, now time.Time, halfLifeDays float64) float64 {
@@ -228,79 +192,6 @@ func weightedResidualsPowerLaw(params, weight, reps []float64, dates []time.Time
 	return sum
 }
 
-// fitLinear fits w = a + b*reps
-func fitLinear(weight, reps []float64, dates []time.Time, now time.Time, halfLifeDays float64) (Estimator, float64, error) {
-	params := []float64{weight[0], 0.0}
-	problem := optimize.Problem{
-		Func: func(x []float64) float64 {
-			return weightedResidualsLinear(x, weight, reps, dates, now, halfLifeDays)
-		},
-	}
-	result, err := optimize.Minimize(problem, params, nil, nil)
-	if err != nil {
-		return nil, 0, err
-	}
-	residual := weightedResidualsLinear(result.X, weight, reps, dates, now, halfLifeDays)
-	return &LinearEstimator{
-		a:          result.X[0],
-		b:          result.X[1],
-		halfLife:   halfLifeDays,
-		modelType:  ModelLinear,
-		residualSum: residual,
-	}, residual, nil
-}
-
-func weightedResidualsLinear(params, weight, reps []float64, dates []time.Time, now time.Time, halfLifeDays float64) float64 {
-	a, b := params[0], params[1]
-	var sum float64
-	for i := range weight {
-		daysAgo := now.Sub(dates[i]).Hours() / 24
-		weightDecay := math.Exp(-math.Ln2 * daysAgo / halfLifeDays)
-		predicted := a + b*reps[i]
-		residual := weight[i] - predicted
-		sum += weightDecay * residual * residual
-	}
-	return sum
-}
-
-// fitExponential fits w = a * exp(b*reps)
-func fitExponential(weight, reps []float64, dates []time.Time, now time.Time, halfLifeDays float64) (Estimator, float64, error) {
-	params := []float64{max(weight), -0.01}
-	problem := optimize.Problem{
-		Func: func(x []float64) float64 {
-			return weightedResidualsExponential(x, weight, reps, dates, now, halfLifeDays)
-		},
-	}
-	result, err := optimize.Minimize(problem, params, nil, nil)
-	if err != nil {
-		return nil, 0, err
-	}
-	residual := weightedResidualsExponential(result.X, weight, reps, dates, now, halfLifeDays)
-	return &ExponentialEstimator{
-		a:          result.X[0],
-		b:          result.X[1],
-		halfLife:   halfLifeDays,
-		modelType:  ModelExponential,
-		residualSum: residual,
-	}, residual, nil
-}
-
-func weightedResidualsExponential(params, weight, reps []float64, dates []time.Time, now time.Time, halfLifeDays float64) float64 {
-	a, b := params[0], params[1]
-	var sum float64
-	for i := range weight {
-		daysAgo := now.Sub(dates[i]).Hours() / 24
-		weightDecay := math.Exp(-math.Ln2 * daysAgo / halfLifeDays)
-		predicted := a * math.Exp(b*reps[i])
-		residual := weight[i] - predicted
-		sum += weightDecay * residual * residual
-	}
-	return sum
-}
-
-// --- Utility Functions ---
-
-// max returns the maximum value in a slice
 func max(slice []float64) float64 {
 	m := slice[0]
 	for _, v := range slice {
@@ -311,3 +202,32 @@ func max(slice []float64) float64 {
 	return m
 }
 
+func extractA(data []timePoint) []float64 {
+	out := make([]float64, len(data))
+	for i, d := range data {
+		out[i] = d.a
+	}
+	return out
+}
+
+func extractB(data []timePoint) []float64 {
+	out := make([]float64, len(data))
+	for i, d := range data {
+		out[i] = d.b
+	}
+	return out
+}
+
+// exponentialSmoothing applies exponential smoothing to a time series.
+func exponentialSmoothing(series []float64, alpha float64) []float64 {
+	if len(series) == 0 {
+		return nil
+	}
+	smoothed := make([]float64, len(series))
+	smoothed[0] = series[0]
+	for i := 1; i < len(series); i++ {
+		smoothed[i] = alpha*series[i] + (1-alpha)*smoothed[i-1]
+	}
+	return smoothed
+}
+